Article pubs.acs.org/IECR
Effect of Shear-Thinning Behavior on Heat Transfer from a Heated Sphere in Yield-Stress Fluids N. Nirmalkar,† R. P. Chhabra,†,* and R. J. Poole‡ †
Department of Chemical Engineering, Indian Institute of Technology, Kanpur, India 208016 School of Engineering, University of Liverpool, Liverpool, L69 3GH, United Kingdom
‡
ABSTRACT: The present study deals with the prediction of drag and forced convection heat transfer behavior of a heated sphere in shear-thinning yield-stress fluids over wide ranges of conditions: plastic Reynolds number, 1 ≤ Re ≤ 100; Prandtl number, 1 ≤ Pr ≤ 100; Bingham number, 10−3 ≤ Bn ≤ 10; and shear-thinning index, 0.2 ≤ n ≤ 1. The momentum and energy equations have been solved numerically together with the Papanastasiou regularization method for viscosity to circumvent the discontinuity inherent in the Herschel−Bulkley constitutive equation. Extensive results on the flow and heat transfer characteristics are presented in order to delineate the influence of the aforementioned dimensionless parameters. Thus, for instance, the flow characteristics are presented in terms of the streamlines, morphology of the yielded/unyielded regions, recirculation length, shear rate magnitude over the surface of the sphere, and drag coefficient. Similarly, heat transfer characteristics are examined in terms of isotherm contours in the close proximity of the sphere and the average Nusselt number as a function of the relevant dimensionless groups. Furthermore, the present results are compared with the available experimental and numerical results in order to establish the reliability and precision of the numerical solution methodology employed in this work. Finally, the average Nusselt number and drag values are correlated in terms of the shear-thinning index (n) and the modified Reynolds number (Re*) via simple expressions, thereby enabling their interpolation for intermediate values of the modified Reynolds number. All else being equal, in addition to Bingham number, shear-thinning behavior of yield stress fluids enhances the rate of heat transfer over and above that observed in Newtonian fluids.
1. INTRODUCTION Yield stress fluids occur frequently in the formulation and processing of numerous engineering products, and examples include foodstuffs such as fruit purees and pastes, personal and health care products like sunscreens, toothpastes, hydrocolloids, emulsions, and foams. Additional examples are found in the context of coating colors and suspensions. Often processing of such fluids entails heat transfer, e.g., during the curing, sterilization, or cooking of foodstuffs.1,2 Furthermore, similar nonisothermal effects also occur during the flow of yield stress fluids in porous media such as during the use of sand pack filters and in the enhanced oil recovery via the use of emulsions and foams. Hence, as such it is important to understand the momentum and heat transfer processes in yield stress fluids in such model configurations. This, in turn, can be used as a basis to model transport processes in complex multiparticle systems including fixed and fluidized bed reactors, porous media flows, etc.3 Besides, such model configurations like momentum and heat transfer from a single sphere also contribute to the overall growth of this subdomain of non-Newtonian fluid mechanics. In a recent study,4 the momentum and forced convection heat transfer characteristics of a heated sphere immersed in Bingham plastic fluids have been studied numerically in the axisymmetric flow regime. Extensive numerical results were reported not only on drag coefficient and Nusselt number but also on the detailed structures of the velocity and temperature fields as functions of the Reynolds number (1 ≤ Re ≤ 100), Prandtl number (1 ≤ Pr ≤ 100), and Bingham number (0 ≤ Bn ≤ 104). While the fluid-like or yielded zones diminish with the increasing value of the Bingham number, the flow remains attached to the surface of the © 2013 American Chemical Society
Figure 1. (a) Schematic of the flow. (b) Computational domain.
sphere even up to Re = 100 at as small a value of the Bingham number as Bn = 1. This is in stark contrast to the case of Newtonian fluids wherein flow detachment occurs at Re ≅ 20−22.5 Therefore, in such visco-plastic fluids, the yield stress acts to Received: Revised: Accepted: Published: 13490
July 4, 2013 August 21, 2013 August 21, 2013 August 21, 2013 dx.doi.org/10.1021/ie402109k | Ind. Eng. Chem. Res. 2013, 52, 13490−13504
Industrial & Engineering Chemistry Research
Article
Figure 2. Representative streamline contours at Re = 1 (top four rows) and Re = 50 (bottom four rows).
(creeping flow) regime. Furthermore, depending upon the density of the sphere (ρs) with respect to that of the fluid (ρ), Tabuteau et al.16 observed three settling regimes. For ρs ≫ ρ, the sphere attained its terminal velocity akin to that in Newtonian fluids. Below a critical value of the sphere density, the sphere came to a complete halt, whereas for the intermediate values of ρs, the sphere never attained a constant terminal velocity. They were able to explain these results with the help of creep tests performed on their test fluids together with some aging effects. Also, their results in the so-called steady-flow regime were in line with previous studies.10,17 No analogous heat transfer results are available for a heated sphere submerged in Herschel−Bulkley fluids. The objectives of the present study are thus 2-fold: First, we examine the combined effects of shear-thinning and fluid
stabilize the flow, whereas the inertial effects tend to destabilize it. Furthermore, the numerically predicted4 values of drag coefficient were found to be in fair agreement with the scant experimental results available in the literature. It is readily acknowledged that most visco-plastic fluids also exhibit varying levels of shear-thinning once the prevailing stress levels exceed the fluid yield stress.6−9 Therefore the flow behavior of the so-called yieldpseudoplastic fluids is frequently approximated much more closely by the familiar Herschel−Bulkley fluid model than a simple Bingham fluid model, e.g., see refs 10−17. While limited numerical and experimental results of drag on spheres10,14−16 in Herschel−Bulkley fluids are available in the literature, most of these are concerned with the development of the stability criterion and drag expressions in the low Reynolds number 13491
dx.doi.org/10.1021/ie402109k | Ind. Eng. Chem. Res. 2013, 52, 13490−13504
Industrial & Engineering Chemistry Research
Article
Figure 3. Representative streamline contours at Re = 70 (top four rows) and Re = 100 (bottom four rows).
inertia on drag and Nusselt number behavior of an isothermal heated sphere. Second, in shear-thinning fluids without yield stress, all else being equal, shear-thinning behavior promotes the rate of convective heat/mass transport over and above that observed in Newtonian fluids for a sphere in the forced-, free-, and mixed-convection regimes.18−22 Intuitively, it appears that while the presence of a yield stress suppresses fluid deformation, shear-thinning viscosity is likely to enhance the extent of convection due to the lowering of the effective fluid viscosity. It is thus useful to elucidate the combined effects of fluid inertia and shear-thinning on the drag and Nusselt number behavior of a sphere in Herschel−Bulkley fluids. Since the bulk of the pertinent literature has been reviewed thoroughly in refs 4 and 8,
it is not repeated here. Suffice it to add here that most of the available drag results are restricted to the so-called creeping-flow regime, and there are no prior results on convective heat transfer for a sphere in Herschel−Bulkley model fluids. This work endeavors to fill this gap in the literature. In particular, the governing partial differential equations (continuity, momentum, and energy) have been solved numerically for the following ranges of conditions: plastic Reynolds number, 1 ≤ Re ≤ 100; Prandtl number, 1 ≤ Pr ≤ 100; Bingham number, 10−3 ≤ Bn ≤ 10; and shear-thinning index, 0.2 ≤ n ≤ 1 (the definitions for these nondimensional parameters are provided in the following section). This work endeavors to elucidate 13492
dx.doi.org/10.1021/ie402109k | Ind. Eng. Chem. Res. 2013, 52, 13490−13504
Industrial & Engineering Chemistry Research
Article
2. PROBLEM FORMULATION Since we are considering the same configuration (Figure 1) as that studied in our recent paper,4 the discussion of the problem statement and governing equations, boundary conditions, etc. is not repeated here. For incompressible fluids, the deviatoric stress tensor can be written as follows: τ = ηγ ̇
(1)
where γ̇ is rate-of-strain tensor. The magnitudes of the rate of deformation tensor and the deviatoric stress tensor, respectively, are given by
Figure 4. Variation of recirculation length (Lr) with shear-thinning index (n).
|γ |̇ =
1 tr(γ 2̇ ) 2
(2)
|τ | =
1 tr(τ 2) 2
(3)
The combined effects of the yield stress and shear-thinning behavior are approximated here by the familiar Herschel− Bulkley fluid model which, in steady simple shear flow, is written in nondimensional form as
the influence of each of these parameters on the momentum and heat transfer characteristics of a heated sphere.
Figure 5. Effect of Reynolds number, Bingham number, and power-law index on morphology of yielded/unyielded regions. 13493
dx.doi.org/10.1021/ie402109k | Ind. Eng. Chem. Res. 2013, 52, 13490−13504
Industrial & Engineering Chemistry Research
Article
nondimensionalized as (T′ − T0)/(Tw − T0). Dimensional considerations suggest that the flow field and drag behavior of the sphere are influenced by three dimensionless parameters, namely, power-law index (n) and Reynolds (Re) and Bingham numbers (Bn), whereas the temperature field and Nusselt number (Nu) show additional dependence on the Prandtl number (Pr). These are defined as follows: Bingham number. Bingham number: Bn =
τY K (V0/d)n
(7)
In the limit of Bn→0, eq 4 reduces to power-law fluid behavior (viscous limit), and on the other hand, in the limit of Bn→∞, it approaches the condition of purely plastic flow. The form of the Bingham number given by eq 7 is also sometimes known as the Oldroyd number in the literature, e.g., see references 11−13 and 29. Reynolds number. For Herschel−Bulkley fluids, it can be defined as follows: ρd nV02 − n K
Re =
(8)
Prandtl number. For Herschel−Bulkley fluids, it is usually expressed as follows: Pr =
(4)
γ ̇ = 0, if |τ | ≤ Bn
(5) 23−27
Much has been written in the literature about the various regularization schemes formulated to treat the discontinuity inherent in the Herschel−Bulkley fluid model; suffice it to add here that the so-called Papanastasiou23 approach is used in the current simulations, following its success in our previous studies.4,28 Within this framework, the Herschel−Bulkley model is rewritten as ⎛ Bn[1 − exp(−m|γ |̇ )] ⎞ τ = ⎜|γ |̇ n − 1 + ⎟γ ̇ |γ |̇ ⎝ ⎠
(9)
These dimensionless groups are obtained by using the representative viscosity given as η ∼ K(V0/d)n−1, thereby disregarding the role of the yield stress in this approach. However, in some circumstances, it is perhaps more appropriate to use the effective viscosity given by ηeff ∼ (τY(d/V0)) + K(V0/d)n−1, which still uses the representative shear rate given by γ̇ ∼ V0/d, but incorporates the contribution of yield stress. Therefore, the resulting modified definitions of the Reynolds number (Re*) and Prandtl number (Pr*) can be written as follows:
Figure 6. Comparison of yield surface predicted by the Papanastasiou regularization (m = 106) and biviscosity model (μy/η ≥ 104) at Re = 100, Bn = 1, and (a) n = 0.5 or (b) n = 0.3.
⎛ Bn ⎞ τ = ⎜|γ |̇ n − 1 + ⎟γ ̇, if |τ | > Bn |γ |̇ ⎠ ⎝
n−1 CK ⎛ V0 ⎞ ⎜ ⎟ k ⎝d⎠
Re* =
Re Bn + 1
Pr * = Pr(Bn + 1)
(10)
In our previous study,4 the governing differential equations subject to the usual boundary conditions (no-slip and constant temperature on the surface of the sphere; uniform flow condition at the inlet plane; slip boundary condition at the confining walls; and zero gradients of velocity and temperature in the axial direction at the outlet) were solved numerically, as described recently.4 Suffice it to say here that the computational mesh and domain employed previously4 were found to be adequate for the results to be free from any numerical artifacts. Finally, the value of m = 106 in eq 3 was found to be satisfactory over the range of conditions spanned here, in line with our previous studies.4,28 Once the flow domain is mapped in terms of the velocity and temperature, one can deduce the values of drag coefficient (CD), recirculation length (Lr), and surface average Nusselt number (Nu) as functions of the Reynolds number (Re), Bingham number (Bn), Prandtl number (Pr), and power-law index (n). This work endeavors to understand these functional relationships. Since extensive validation of the numerical solution methodology and of the choice of parameters (Lu = 20, Ld = 30, H = 30, m = 106) has been presented in ref 4, no additional benchmark comparisons are included here.
(6)
In eq 6, the exponent, m (nondimensionalized using d/V0), controls the stress growth. For stress levels below the fluid yield stress, it allows a finite value of stress corresponding to vanishingly small shear rates, whereas beyond the yield stress, the applied stress varies with the shear rate similar to the power-law type behavior. On the other hand, when m→0, it would correctly approach the limit of power-law fluids (Bn = 0), and similarly, in the limit of m→∞, it leads to the true Herschel−Bulkley model behavior. The aforementioned equations are rendered dimensionless using sphere diameter, d; free stream fluid velocity, V0; K(V0/d)n; d/V0; and V0/d as scaling variables for length, velocity, stress, time, and shear rate respectively. The temperature was 13494
dx.doi.org/10.1021/ie402109k | Ind. Eng. Chem. Res. 2013, 52, 13490−13504
Industrial & Engineering Chemistry Research
Article
Figure 7. Distribution of shear rate magnitude over the sphere surface at Re = 1.
3. RESULTS AND DISCUSSION In this study, the drag and forced convection heat transfer behavior of a heated sphere in shear-thinning yield stress fluids modeled using the Herschel−Bulkley equation have been examined in detail over wide ranges of conditions: 1 ≤ Re ≤ 100; 1 ≤ Pr ≤ 100; 10−3 ≤ Bn ≤ 10; and 0.2 ≤ n ≤ 1. Extensive new results regarding streamline patterns, recirculation length, morphology of yielded/unyielded regions, kinematics of flow, drag coefficient, isotherm contours, and average Nusselt number are presented and discussed spanning the aforementioned range of the pertinent dimensionless parameters. While no definitive information is available regarding the flow regimes, it is assumed to be axisymmetric here as the maximum value of the Reynolds number is only 100 in this study, which is significantly below the value of Re = 250 reported for Newtonian fluids.5 Besides, the yield stress also acts to stabilize the flow. Thus, on either count, the assumption of the axisymmetric flow is well justified here. 3.1. Streamline Contours. Representative streamline contours are shown in Figures 2 and 3 for a range of values of Reynolds number, Bingham number, and shear-thinning index. At Re = 1, and irrespective of the values of the power-law index (n) and of the Bingham number (Bn), no flow separation is observed. Such flow attachment is simply due to the negligible fluid inertia at such low Reynolds numbers. On the other hand, a visible recirculation zone is seen to have formed at Re = 50 for values of Bingham number Bn ≤ 1 irrespective of the value of shear-thinning index, n. For a fixed value of the Reynolds number, the recirculation region is seen to progressively diminish, eventually disappearing altogether with the increasing value of the Bingham number; e.g., see the results for Bn = 1 at Re = 50
where a very small wake is present which disappears completely at Bn = 10, similar to that seen in Bingham plastic liquids.4 The effect of the power-law index is seen to be rather negligible under these conditions. Intuitively, it appears that with the increasing Reynolds number, higher values of yield stress, i.e., Bingham number, would be needed to suppress the propensity for flow separation. This conjecture is well borne out by the results shown in Figure 3 at Re = 70 and Re = 100. This finding is in line with that reported by Mossaz et al.29 for the flow past a circular cylinder. Furthermore, for Bingham numbers greater than unity, i.e., Bn > 1, no flow separation was seen to occur for any value of the shear-thinning index, n, spanned here even at Re = 100, as shown in Figure 4. Naturally, this is due to the yield stress of the fluid, which tends to suppress fluid deformation and acts to stabilize the flow field. The size of the wake is described here in terms of the recirculation length (Lr; measured from the rear stagnation point of the sphere) in Figure 4. Inspection of this figure shows the following trends. First, in the limit of negligible yield stress, the recirculation length in shear-thinning fluids (n < 1) is seen to be smaller than that in Newtonian fluids (n = 1) under most conditions. On the other hand, in visco-plastic fluids, the recirculation length (Lr) shows a reverse dependence on the shearthinning index (n). We believe that this switch over is due to the interaction between the two nonlinear forces (inertial and viscous), both of which scale differently with velocity. For instance, the rate of fluid deformation is at a maximum close to the surface of the sphere and it decreases rapidly due to both the yield stress and shear-thinning behavior. Both tend to suppress the degree of recirculation in the separated flow region. 13495
dx.doi.org/10.1021/ie402109k | Ind. Eng. Chem. Res. 2013, 52, 13490−13504
Industrial & Engineering Chemistry Research
Article
Figure 8. Distribution of shear rate magnitude over the sphere surface at Re = 100.
3.2. Morphology of Yielded/Unyielded Zones. Representative results on the yielded/unyielded zones around the sphere are presented in Figure 5 for a range of values of the Reynolds number, Bingham number, and shear-thinning index. In this study, for a fixed value of the Bingham number, the fluidlike zone is seen to expand on increasing the Reynolds number irrespective of the value of shear-thinning index. This increase is obviously due to the increasing importance of the inertial forces and the diminishing viscous forces. Similarly, for a fixed value of the Reynolds number, the size of the fluid-like region shrinks, and it is confined to the proximity of the sphere upon increasing the Bingham number. Furthermore, the size of the fluid-like zone in Herschel−Bulkley fluids is seen to be somewhat smaller than that in Bingham plastic fluids except at Re = 1 and 10. This is probably due to the rapid spatial decay of the velocity field in shearthinning fluids. The so-called “polar cap” at the front stagnation point disappears altogether at finite values of the Reynolds number due to inertial effects, whereas the one in the rear of the sphere is seen to grow in size with the Reynolds number. Intuitively, as the yielded zones shrink with the increasing Bingham number, it is thus fair to postulate that the influence of shear-thinning index will progressively vanish with the increasing Bingham number. Indeed, this behavior was seen to occur at Bn = 10, where the whole body of fluid is seen to be virtually unyielded. In summary, fluid-like regions shrink with the increasing Bingham number. On the other hand, fluid-like regions increase in their spatial extent with increasing Reynolds number. Also, for shear-thinning yield-stress fluids, the fluid-like regions are seen to be somewhat smaller than that in Bingham plastic fluids otherwise under identical conditions: this effect is
especially noticeable at Bn = 1 and Re = 50 and 100. Finally, included in Figure 6 are also the predictions of the widely used biviscous model approach (shown as broken lines and for μy/η ≥ 104), and the two predictions are seen to be fairly close to each other, thereby inspiring confidence in the reliability of the new numerical results reported herein. 3.3. Flow Kinematics. The kinematics of the flow are analyzed here in terms of the distribution of the magnitude of shear rate along the surface of the sphere (Figures 7 and 8). The maximum rate of deformation in the fluid is seen to occur at the equator, i.e., θ = 90° for Re = 1 as shown in Figure 7. On the other hand, at Re = 100, this maxima is seen to move toward the front stagnation point and occurs in between θ = 0° and θ = 90°. For a fixed value of the Reynolds number, the peak value of the shear rate increases with Bingham number for all values of shearthinning index. This is simply due to the fact that the yielded zone is confined to a very thin layer, and due to the imposition of the no-slip conditions on the surface of the sphere, gradients become very steep. Similarly, all else being equal, shear-thinning fluids also experience a somewhat larger rate of deformation in the fluid, and this can safely be attributed to the further thinning of the momentum boundary layer on account of shear-thinning behavior. Some further insights are provided by investigating the velocity and shear rate contour plots (Figure 9a,b). As noted earlier, owing to the presence of the yield stress, fluid-like zones closely follow the surface of the sphere, and therefore a relatively high magnitude of velocity is seen to occur in the vicinity of the sphere, as shown in Figure 9a. Moreover, shear-thinning behavior is also seen to suppress the flow field in addition to yield stress. Similar trends are observed in terms of the shear rate magnitude 13496
dx.doi.org/10.1021/ie402109k | Ind. Eng. Chem. Res. 2013, 52, 13490−13504
Industrial & Engineering Chemistry Research
Article
Figure 9. Contour diagrams: (a) velocity magnitude (top two rows) and (b) shear rate magnitude (bottom two rows).
another factor, α, to obtain improved correlation of their data as follows:
contours as shown in Figure 9b. In a nutshell, while the influence of the shear-thinning fluid behavior and yield-stress go handin-hand in suppressing the yielded (or fluid-like) domains, this tendency is opposed by the increasing inertial forces, i.e., increasing Reynolds number. However, as the yielded regions shrink in size, the velocity and temperature gradients sharpen, which, in turn, should lead to varying levels of augmentation both in the drag and Nusselt number of the sphere. 3.4. Drag Coefficients. As noted earlier, the drag coefficient represents the net hydrodynamic force exerted by the fluid on the sphere in the direction of flow. In this work, the drag coefficient is presented as a function of the modified Reynolds number and shear-thinning index. It is convenient to treat the low (creeping flow) and high Reynolds number results separately. 3.4.1. Creeping-Flow Region. At low Reynolds numbers, the use of the modified Reynolds number (Re*) leads to the reconciliation of the experimental results for different values of the power-law index.10,16 However, it needs to be emphasized here that both Atapattu et al.10 and Tabuteau et al.16 introduced
Re* =
Re 1 + αBn
(11)
Naturally, in the limit of α = 1, eq 11 coincides with the definition of Re* given by eq 10. Using a slip-line analysis, Ansley and Smith30 proposed a value of α = (7π/24) which is not only close to unity but is also within ±10% of its numerical value.13 While Atapattu et al.10 were able to correlate their experimental results (0.43 ≤ n ≤ 0.84) adequately using α = 1, Tabuteau et al.16 argued that the theoretical value of α = 0.823 was more appropriate, albeit their results were limited to a single value of n = 0.5. Furthermore, Atapattu et al.7 incorporated explicitly a correction factor X(n) for the effect of power-law index n, which is available in the literature based on several numerical studies, e.g., see ref 8. Thus, Atapattu et al.10 proposed the following expression for drag in the creeping-flow region: 13497
dx.doi.org/10.1021/ie402109k | Ind. Eng. Chem. Res. 2013, 52, 13490−13504
Industrial & Engineering Chemistry Research
Article
Figure 12. Effect of modified Reynolds number (Re*) on drag coefficient (CD) at finite Reynolds numbers.
Figure 10. Effect of modified Reynolds number (Re*) on drag coefficient (CD) in the limit of creeping flow.
Table 2. Values of Constants in eq 13
Table 1. Range of Shear-Thinning Index, n, in Experimental Studies ref
n 0.43−0.84 0.5 1 0.55 0.56−0.85 0.48−0.5
CD =
24X(n) Re*
a
b
c
29.544 30.960 32.568 34.248 35.952 37.680 39.454 41.237 43.032
0.027 0.022 0.016 0.012 0.009 0.007 0.005 0.004 0.003
0.724 0.806 0.907 1.008 1.098 1.177 1.247 1.310 1.368
plotted in accordance with eq 12. It is clearly seen that this form of presentation consolidates the present drag results for all values of n and up to about Re* = 1; the value of CDRe* is nearly constant and is within ±28.5% of 24X(n). This level of deviation is in line with the ±25% uncertainty of such experiments reported by Atapattu et al.10 Figure 11 shows a comparison between the present numerical predictions and the scant experimental results available in the literature where the correspondence is seen to be fair. The range of shear-thinning index n embraced in the experimental studies is summarized in Table 1. For the sake of comparison, the creeping flow numerical predictions14 are also included in this figure, which are seen to be in excellent agreement with the present low Reynolds number results. 3.4.2. Finite Reynolds Number Region. Figure 12 shows the present drag results plotted as a function of the modified Reynolds number, Re*, and it is clearly seen that this form of representation leads to a family of curves depending upon the value of the shear-thinning index (n). Broadly, the larger is the value of n, the higher is the drag beyond the creeping-flow regime. Conversely, drag reduction occurs in Herschel−Bulkley fluids with reference to that in Bingham plastic fluids. In the creeping flow regime, the pressure forces are in equilibrium with the viscous forces (the latter also includes the effect of the fluid yield stress). Beyond the creeping-flow regime, inertial forces also exert varying levels of influence in determining the net value of drag exerted on the sphere, especially as the inertial and viscous forces scale differently with velocity. Similar problems were also encountered by Valentik and Whitmore32 and Pazwash
Figure 11. Comparison of the present drag results with the literature experimental results in the creeping flow region.
Atapattu et al.10 Tabuteau et al.16 Chafe and de Bruyn31 Hariharaputhiran et al.17 Sen34 Beaulne and Mitsoulis14 (numerical results)
n 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
(12)
Though the form of eq 12 is difficult to justify on a sound physical basis, it was found to be satisfactory (±25%) to correlate their results spanning ranges of conditions as 0.43 ≤ n ≤ 0.84, 0 ≤ Bn ≤ 280, and Re ≤ 0.36. Figure 10 shows the present results 13498
dx.doi.org/10.1021/ie402109k | Ind. Eng. Chem. Res. 2013, 52, 13490−13504
Industrial & Engineering Chemistry Research
Article
Figure 13. Representative isotherm contours at Re = 1 (top four rows) and Re = 50 (bottom four rows).
While the values of the fitted constants in eq 13 are given in Table 2, the overall resulting mean and maximum deviations are 4% and 6.8%, respectively, without any discernible trends. 3.5. Isotherm Contours. Typical isotherm contours in the vicinity of the sphere are shown in Figures 13 and 14 for representative values of the Reynolds number, Bingham number, Prandtl number, and power-law index. For a fixed value of the Prandtl number and at low Reynolds numbers when convection
and Robertson33 in correlating their high Reynolds number sphere drag results in Bingham plastic fluids. Finally, before leaving this section, the present numerical results are reasonably well approximated by the following expression (shown by lines in Figure 12): CD =
a (1 + bRe*c ) Re*
(13) 13499
dx.doi.org/10.1021/ie402109k | Ind. Eng. Chem. Res. 2013, 52, 13490−13504
Industrial & Engineering Chemistry Research
Article
Figure 14. Representative isotherm contours at Re = 70 (top four rows) and Re = 100 (bottom four rows).
is weak, heat transfer is dominated by conduction. Under these conditions, the isotherms are seen to be more or less similar in the fore and aft regions, and this trend persists in highly shearthinning fluids and/or at high Bingham numbers also due to the vanishing extent of yielded zones. On the other hand, all else being equal, at high Reynolds numbers, isotherms are seen to cluster in the vicinity of the sphere especially in the front of the sphere and are increasingly elongated in the rear of the sphere
due to the increasing convection effects. Similar trends are observed with the increasing Bingham and Prandtl numbers. This distortion is obviously due to the thinning of the thermal and hydrodynamic boundary layers at a high Reynolds number, Prandtl number, and Bingham number. Furthermore, the clustering of isotherms is seen to be somewhat more severe in Herschel−Bulkley fluids than that in Bingham plastic fluids (n = 1), thereby suggesting some enhancement in heat transfer on account 13500
dx.doi.org/10.1021/ie402109k | Ind. Eng. Chem. Res. 2013, 52, 13490−13504
Industrial & Engineering Chemistry Research
Article
Figure 15. Effect of shear-thinning index (n) on the average Nusselt number at Re = 100.
Table 3. Values of a1 in eq 15 n
a1
average % error
maximum % error
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
13.35 10.22 7.870 6.070 4.701 3.977 3.122 2.452 2.290
11.50 12.27 12.55 12.95 12.97 12.88 13.50 13.97 15.76
26.03 28.01 28.45 27.50 28.98 28.94 26.92 35.60 40.54
in Figure 15), Prandtl number, and Bingham number. Naturally, this functional relationship Nu = f(Re, Pr, Bn, n) will result in a family of curves depending upon the values of these parameters. Following the success of the modified Reynolds (Re*) and Prandtl (Pr*) number coordinates by incorporating the influence of the Bingham number, this functional relationship can be reduced to Nu = f (Re*, Pr*, n). This format will still yield a family of curves corresponding to each value of Pr* and n. Further consolidation of the results is achieved here by employing the familiar Colburn factor, jH, as shown in our previous studies.4,22 The Colburn factor is defined as follows:
Figure 16. Dependence of jH factor on the modified Reynolds number (Re*).
of shear-thinning behavior over and above that seen in Bingham plastic fluids. 3.6. Average Nusselt Number. The influence of the shear-thinning index on the average Nusselt number is shown in Figure 15 at Re = 100. As expected, it shows a monotonic increase with the decreasing value of shear-thinning index (n). Furthermore, as postulated here, the average Nusselt number shows a positive dependence on the Reynolds number (not shown
jH =
Nu = f (Re*, n) Re*Pr *1/3
(14)
Intuitively, it appears that this format will still yield a family of curves for each value of n, as is borne out by Figure 16 where the complete set of the present numerical results is plotted in accordance with eq 14. On the other hand, as expected, the use 13501
dx.doi.org/10.1021/ie402109k | Ind. Eng. Chem. Res. 2013, 52, 13490−13504
Industrial & Engineering Chemistry Research
Article
Table 4. Comparison between the Present Results for Bn = 10−3 and the Power-Law Model Predictions power-law fluids Song et al.20,36 n
Re
Pr
Nuavg
0.2
100
1 5 10 1 5 50 100 1 5 10 20 1 5 10 20 1 5 50 100 1 5 10 20 1 5 10 20 1 5 10 100 1 5 10 20
10.01 19.05 25.25 3.155 5.735 9.307 11.511 6.559 10.95 13.61 16.711 8.611 14.726 18.231 21.984 3.078 5.379 8.53 10.471 6.227 10.112 12.458 15.211 8.1 13.477 16.569 19.927 3.0841 4.4120 5.1230 9.6750 7.8970 12.523 15.234 18.354
0.6
5
50
100
0.8
5
50
100
1
5
100
Herschel−Bulkley fluids Dhole et al.18,35
CD
Nuavg
present CD
0.5321
7.306
1.219
0.795
7.23
1.402
0.951
7.143
1.023
3.139 5.667 9.063 11.183 6.512 10.778 13.374 16.482 8.556 14.652 18.46 22.62 3.081 5.389 8.554 10.551 6.234 10.14 12.556 15.59 8.11 13.567 16.97 20.65 2.9243 4.3460 5.2359 10.216 7.7180 12.543 15.563 19.369
6.810
1.170
0.760
7.040
1.370
0.920
7.080
1.060
Nuavg 10.02 19.16 25.28 3.156 5.737 9.320 11.589 6.559 10.966 13.716 17.174 8.614 14.853 18.770 23.699 3.079 5.381 8.541 10.536 6.227 10.125 12.544 15.590 8.101 13.575 17.001 21.311 3.0865 4.4321 5.2351 9.6876 7.9098 12.543 15.243 18.378
CD 0.5255
7.345
1.223
0.789
7.265
1.4045
0.9402
7.126
1.022
advanced by Song et al.36 and thus are not repeated here. Intuitively, it appears that for given values of the power-law index and Reynolds number, there will be a critical Bingham number above which the effects of the yield stress begin to influence the results. At this juncture, it is important to make two observations here regarding the results reported herein. The use of the regularized Herschel−Bulkley model is really tantamount to using the generalized Newtonian fluid (GNF) approach, which hinges on the existence of a scalar viscosity function and as such does not entertain the existence of the yield stress. The fact that the two regularization schemes, namely, biviscosity and exponential forms of the Herschel−Bulkley and/or Bingham plastic model, lead to virtually identical results not only inspires confidence in the reliability and accuracy of the results reported herein but also suggests that both of these regularization schemes can be used to obviate the discontinuous flow behavior exhibited by yield stress fluids. On the other hand, this approach does not rule out the possibility of obtaining similar results using another appropriate GNF model with the matching viscosity and rate of change of viscosity (dη/dγ̇) as that used here. Second, the thermo-physical properties (i.e., heat capacity, C; thermal conductivity, k; density,
of jH does consolidate results for scores of values of Pr*. The dependence of the power-law index is best approximated by the expression a1 jH = (15) Re*2/3 The values of a1(n) in eq 15 together with the resulting average and maximum percentage errors are reported in Table 3. The effect of the power-law index is seen to be quite dramatic here, and it is thus possible to realize significant enhancement in the value of mean Nusselt number in highly shear-thinning fluids even in yield stress fluids. Finally, before leaving this section, it is worthwhile to examine the present results in the limit of Bn→0. Intuitively, it appears that, in this limit, the present values of drag coefficient and Nusselt number should approach the corresponding values in power-law fluids. Table 4 shows such a comparison between the present results for Bn = 10−3 and that in power-law fluids. Clearly, the close correspondence between the present results and that of Song et al.19,36 indicates little or no effect of such a low value of the Bingham number. The reasons for the deviations on the order of 10% in drag values with that of Dhole et al.35 have been 13502
dx.doi.org/10.1021/ie402109k | Ind. Eng. Chem. Res. 2013, 52, 13490−13504
Industrial & Engineering Chemistry Research
Article
ρ; viscosity, η; or τY, n, K, etc.) are assumed to be temperatureindependent. Among all these properties, perhaps the viscosity is expected to be the most temperature sensitive. While in principle, there is no conceptual difficulty in incorporating this effect into the present numerical analysis, much confusion, however, exists in the literature regarding the effect of temperature on the values of τY, K, and n as noted recently.4 Therefore, more work is needed to delineate the effect of temperature on the Herschel− Bulkley model parameters to develop functional dependence between them and temperature before incorporating these effects into the numerical solutions such as that reported herein.
m = Growth rate parameter, dimensionless n = Shear-thinning index, dimensionless ns = Unit normal vector, dimensionless Nu = Average Nusselt number, dimensionless Pr = Prandtl number, dimensionless Pr* = Modified Prandtl number (≡ Pr(1 + Bn)), dimensionless Re = Reynolds number, dimensionless Re* = Modified Reynolds number (≡ Re/(1 + Bn)), dimensionless T = Temperature of fluid (≡(T′ − T0)/(Tw − T0)), dimensionless T′ = Temperature of fluid, K Tw = Sphere surface temperature, K T0 = Free stream temperature, K V0 = Free stream velocity, m·s−1 r,z = Polar coordinates, dimensionless X(n) = Drag correction factor, dimensionless
4. CONCLUSIONS In this study, the momentum and heat transfer characteristics from a heated sphere submerged in Herschel−Bulkley fluids have been examined numerically in the following ranges of conditions: Reynolds number, 1 ≤ Re ≤ 100; Prandtl number, 1 ≤ Pr ≤ 100; Bingham number, 10−3 ≤ Bn ≤ 10; and shear-thinning index, 0.2 ≤ n ≤ 1. Extensive results regarding streamlines, yielded/ unyielded regions, recirculation length, drag coefficient, isotherms, and Nusselt number are provided as functions of the pertinent dimensionless groups. All else being equal, shearthinning behavior promotes the rate of heat transfer even in the presence of yield stress, which is consistent with that seen in power-law fluids in all regimes of heat transfer. Broadly, both shear-thinning and yield-stress (Bn) act to restrict the size of yielded-zones whereas the fluid inertia (Reynolds number) acts in the opposite manner. The former sharpens the velocity and temperature gradients on the surface of the sphere which augments the values of drag coefficient and Nusselt number. Also, the flow separation is deferred to higher values of the Reynolds number in Herschel−Bulkley fluids than that in fluids without a yield stress otherwise under identical conditions. The present predictions of drag are seen to be in line with the scant experimental results available in the literature. Finally, the numerical drag and Nusselt number values have been correlated using simple forms of expressions, thereby enabling their estimation in a new application.
■
Greek symbols
α = Constant in eq 11, dimensionless γ̇ = Rate of strain tensor, dimensionless η = Reference viscosity, Pa·s ηeff = Effective viscosity, Pa·s θ = Position on the surface of sphere, deg μY = Yielding viscosity, Pa·s ρ = Density of the fluid, kg·m−3 ρs = Density of the sphere, kg·m−3 τ = Extra stress tensor, dimensionless τY = Yield stress, Pa Subscripts
■
0 = Inlet condition w = Sphere surface condition
REFERENCES
(1) Laba, D. Rheological Properties of Cosmetics and Toiletries; MarcelDekker: New York, 1993. (2) Steffe, J. F. Rheological Methods in Food Process Engineering; Freeman Press: East Lansing, MI, 1996. (3) Chhabra, R. P.; Comiti, J.; Machac, I. Flow of Non-Newtonian Fluids in Fixed and Fluidized Beds. Chem. Eng. Sci. 2001, 56, 1. (4) 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. Res. Chem. 2013, 52, 6848. (5) Clift, R.; Grace, J. R.; Weber, M. E. Bubbles, Drops and Particles; Academic Press: New York, 1978. (6) Bird, R. B.; Dai, G. C.; Yarusso, B. J. The Rheology and Flow of Viscoplastic Materials. Rev. Chem. Eng. 1983, 1, 1. (7) Barnes, H. A.; Walters, K. The Yield Stress Myth? Rheol. Acta 1985, 24, 323. (8) Chhabra, R. P. Bubbles, Drops, and Particles in Non-Newtonian Fluids, 2nd ed.; CRC Press: Boca Raton, FL, 2006. (9) Chhabra, R. P.; Richardson, J. F. Non-Newtonian Flow and Applied Rheology: Engineering applications, 2nd ed.; Butterworth-Heinemann: Oxford, U.K., 2008. (10) Atapattu, D. D.; Chhabra, R. P.; Uhlherr, P. H. T. Creeping Sphere Motion in Herschel−Bulkley Fluids: Flow Field and Drag. J. Non-Newtonian Fluid Mech. 1995, 59, 245−265. (11) Jossic, L.; Magnin, A. Drag and Stability of Objects in a Yield Stress Fluid. AIChE J. 2001, 47, 2666. (12) Jossic, L.; Magnin, A. Drag of an Isolated Cylinder and Interactions between Two Cylinders in Yield Stress Fluids. J. NonNewtonian Fluid Mech. 2009, 164, 9. (13) Tokpavi, D. L.; Jay, P.; Magnin, A. Interaction between Two Circular Cylinders in Slow Flow of Bingham Viscoplastic fluid. J. NonNewtonian Fluid Mech. 2009, 157, 175.
AUTHOR INFORMATION
Corresponding Author
*Tel.: 0091 512 2597393. Fax: 0091 512 2590104. E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
NOMENCLATURE a, b, c = Fitting constants in eq 13, dimensionless a1 = Fitting constant in eq 15, dimensionless Bn = Bingham number, dimensionless C = Thermal heat capacity of fluid, J·kg−1·K−1 CD = Drag coefficient ≡ (2FD/ρV02((π/4)d2)), dimensionless d = Diameter of the sphere, m FD = Drag force on the sphere, N h = Heat transfer coefficient, W·m−2·K−1 H = Lateral height of the domain, dimensionless jH = Colburn factor (eq 15), dimensionless k = Thermal conductivity of fluid, W·m−1·K−1 K = Consistency index in Herschel−Bulkley model, Pa·sn Ld = Downstream distance, dimensionless Lr = Recirculation length, dimensionless Lu = Upstream distance, dimensionless 13503
dx.doi.org/10.1021/ie402109k | Ind. Eng. Chem. Res. 2013, 52, 13490−13504
Industrial & Engineering Chemistry Research
Article
(14) Beaulne, M.; Mitsoulis, E. Creeping Motion of a Sphere in TubesFilled with Herschel−Bulkley Fluids. J. Non-Newtonian Fluid Mech. 1997, 72, 55. (15) Putz, A. M.; Burghelea, T. I.; Frigaard, I. A.; Martinez, D. M. Settling of an Isolated Spherical Particle in a Yield Stress Shear-Thinning Fluid. Phys. Fluids 2008, 20, 033102. (16) Tabuteau, H.; Coussot, P.; de Bruyn, J. Drag Force on a Sphere in Steady Motion through a Yield Stress Fluid. J. Rheol. 2007, 51, 125. (17) Hariharaputhiran, M.; Shankar Subramanian, R.; Campbell, G. A.; Chhabra, R. P. The Settling of Spheres in a Viscoplastic Medium. J. NonNewtonian Fluid Mech. 1998, 79, 87. (18) 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. (19) Song, D.; Gupta, R. K.; Chhabra, R. P. Effect of Blockage on Heat Transfer from a Sphere in Power-Law Fluids. Ind. Eng. Res. Chem. 2010, 49, 3849. (20) 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. (21) Prhashanna, A.; Chhabra, R. P. Free Convection in Power-Law Fluids from a Heated Sphere. Chem. Eng. Sci. 2010, 65, 6190. (22) Nirmalkar, N.; Chhabra, R. P. Mixed Convection from a Heated Sphere in Power-Law Fluids. Chem. Eng. Sci. 2013, 89, 49. (23) Papanastasiou, T. C. Flow of Materials with Yield. J. Rheol. 1987, 31, 385. (24) Frigaard, I. A.; Nouar, C. On the Usage of Viscosity Regularisation Methods for Visco-Plastic Fluid Flow Computation. J. Non-Newtonian Fluid Mech. 2005, 127, 1. (25) Liu, B. T.; Muller, S. J.; Denn, M. M. Convergence of a Regularisation Method for Creeping Flow of a Bingham Material about a Rigid Sphere. J. Non-Newtonian Fluid Mech. 2002, 102, 179. (26) Glowinski, R.; Wachs, A. On the Numerical Simulation of Viscoplastic Fluid Flow. In Handbook of Numerical Analysis; Ciarlet, P. G., Ed.; Elsevier: North-Holland, 2011; pp 483−717. (27) O’Donovan, E. J.; Tanner, R. I. Numerical Study of the Bingham Squeeze Film Problem. J. Non-Newtonian Fluid Mech. 1984, 15, 75. (28) Nirmalkar, N.; Chhabra, R. P.; Poole, R. J. On Creeping Flow of a Bingham Plastic Fluid Past a Square Cylinder. J. Non-Newtonian Fluid Mech. 2012, 171−172, 17. (29) Mossaz, S.; Jay, P.; Magnin, A. Criteria for the Appearance of Recirculating and Non-Stationary Regimes behind a Cylinder in a Viscoplastic Fluid. J. Non-Newtonian Fluid Mech. 2010, 165, 1525. (30) Ansley, R. W.; Smith, T. N. Motion of Spherical Particles in a Bingham Plastic. AIChE J. 1965, 13, 1193. (31) Chafe, N. P.; de Bruyn, J. R. Drag and Relaxation in a Bentonite Clay Suspension. J. Non-Newtonian Fluid Mech. 2005, 131, 44. (32) Valentik, L.; Whitmore, R. L. The Terminal Velocity of Spheres in Bingham Plastics. Br. J. Appl. Phys. 1965, 16, 1197. (33) Pazwash, H.; Robertson, J. M. Forces on Bodies in Bingham Fluids. J. Hydraulic Res. 1975, 13, 35. (34) Sen, S. An Experimental Study of the Influence of Yield Stress and Fluid Rheology on Drag Coefficients of Spheres. MS Dissertation, Brigham Young University, Provo, UT, 1984. (35) Dhole, S. D.; Chhabra, R. P.; Eswaran, V. Flow of Power-Law Fluids past a Sphere at Intermediate Reynolds Numbers. Ind. Eng. Chem. Res. 2006, 45, 4773. (36) Song, D.; Gupta, R. K.; Chhabra, R. P. Drag on a Sphere in Poiseuille Flow of Shear-Thinning Power-Law Fluids. Ind. Eng. Res. Chem. 2011, 50, 13105.
13504
dx.doi.org/10.1021/ie402109k | Ind. Eng. Chem. Res. 2013, 52, 13490−13504