Crossflow of viscoplastic materials through tube bundles

J. Non-Newtonian Fluid Mech., 75 (1998) 43 – 54 Crossflow of viscoplastic materials through tube bundles Angela O. Nieckele, Mônica F. Naccache, Paulo R. Souza Mendes * Department of Mechanical Engineering Pontifı́cia Uni6ersidade Católica, Rio de Janeiro, RJ 22453 -900, Brazil Received 1 May 1997 Abstract The flow of viscoplastic materials through staggered arrays of tubes is analyzed. The mechanical behavior of the materials is assumed to obey the generalized Newtonian liquid (GNL) model, with a viscosity function given by the biviscosity law. The governing equations of this flow are solved numerically using a finite-volume method with a non-orthogonal mesh. For a representative range of the relevant parameters, results are presented in the form of velocity, pressure and viscosity fields. The pressure drop is also given as a function of rheological and geometric parameters. © 1998 Elsevier Science B.V. Keywords: Viscoplastic materials; Generalized Newtonian liquid; Finite-volume method; Geometric parameters 1. Introduction Crossflows through tube bundles are found in some common engineering situations, such as the shell-side flow in shell-and-tube heat exchangers. Also, these flows are sometimes employed as a two-dimensional idealization of flows through porous media [1]. In heat exchanger design, important information is the relationship between flow rate and pressure drop. If the fluid flowing externally to the tubes is Newtonian, pressure drop and heat transfer information abounds in the literature, and is reviewed in [2]. For non-Newtonian liquids, some work is found for flow of pseudoplastic power-law liquids around a cylinder (e.g. [3,4]). For flow past tube bundles of viscoelastic polymeric solutions, pressure drop data are available [5,1]. However, no information is available in the open literature for viscoplastic materials, as far as the authors know. Bird et al. [6] presented an overview of the rheology and flow of viscoplastic materials, and analyzed some simple flow situations using a generalized Newtonian liquid (GNL) model with a Bingham plastic viscosity function. A number of authors discussed the difficulty in using the * Corresponding author. 0377-0257/98/$19.00 © 1998 Elsevier Science B.V. All rights reserved. PII S 0 3 7 7 - 0 2 5 7 ( 9 7 ) 0 0 0 7 9 - 7 A.O. Nieckele et al. / J. Non-Newtonian Fluid Mech. 75 (1998) 43–54 44 von Mises yield criterion to analyze complex flows numerically. Essentially two types of modification of the Bingham plastic viscosity function have been proposed to simplify computation, namely, the biviscosity model [7–9], and Papanastasiou’s model [10]. Both modifications have been used successfully in analyses and numerical simulations of different complex flows (e.g. [11 –16]). Moreover, Lipscomb and Denn [7] observed that yielding and flow must occur everywhere in complex flows in confined geometries, which is generally inconsistent with the classic Bingham plastic model. Wilson [14] showed that, for suitably large yield stresses, yield surfaces can exist in confined complex geometries if the biviscosity law is employed, even when the Bingham plastic limit is approached. Piau [16] explained that, if some deformation in the plug-flow region is allowed (either elastic or viscous), yield surfaces are possible whenever there are regions of deformation (or deformation rate) low enough as to require stress levels below the yield stress to be realized. The present research is concerned with low Reynolds number flows of viscoplastic materials past tube bundles in staggered arrangement, as shown in Fig. 1. Pressure drop information is obtained by numerical integration of the governing equations. In addition, velocity, pressure and viscosity fields are presented and discussed. The results reported here are of practical significance for the food, petroleum, paper and other important industries that deal with viscoplastic materials. 2. The analysis In flows of viscoplastic materials, wall slip is often observed, but for rough solid surfaces the no-slip boundary condition is plausible. In this analysis, we assume that the material sticks perfectly to the solid boundaries. Fig. 1 shows the two-dimensional geometry considered. The flowing material is assumed to be incompressible, so that the equation of mass conservation reduces to div u =0 (1) where u is the velocity vector field. Fig. 1. The geometry. A.O. Nieckele et al. / J. Non-Newtonian Fluid Mech. 75 (1998) 43–54 45 The momentum equation is, for steady flow and negligible external forces, r div(uu) = − grad p+ div t (2) where p is the pressure, t the extra-stress and r the mass density. The constitutive equation chosen to relate the extra-stress to the flow kinematics is the GNL model: t= h(g; )g; (3) where h(g; ) is the viscosity function, g; grad u +(grad u)T is the rate-of-deformation tensor, and g; 1/2 tr(g; 2) is a measure of its magnitude. A Bingham plastic is a GNL material that has a viscosity function of the type h= ! t 0/g; +mp, if t \t 0 otherwise (4) where mp is the so called plastic viscosity, t0 the yield stress, and t 1/2 tr(t 2) a measure of the magnitude of t. As mentioned above, the Bingham plastic viscosity function is not simple to model numerically, and may also lead to inconsistencies in complex flows [7]. The usual approach is to replace it with another viscosity function which preserves its main features. Here we choose the biviscosity model [7–9,15], and adopt the parameter values recommended by Beverly and Tanner [13]: h= ! t 0/g; +mp, if g; \ g; c m otherwise (5) where m = 1000mp; g; c = t0 t0 = m − mp 999mp (6) By examining Fig. 1, it can be observed that the geometry repeats itself, i.e. it can be divided into identical cells. It has been noted that after a few rows in the downstream direction, the flow also repeats itself, attaining the so-called condition of periodically fully developed flow [17]. Rather than solving the problem for several rows of rods, exploring this periodicity concept allows the flow domain to be dramatically reduced, as indicated in Fig. 1, yielding much more accurate numerical solutions for a given number of grid points. The periodically fully developed flow condition implies that the flow field at the exit of the domain is exactly the same as at the inlet, and that the pressure distribution can be decomposed as p(x, y)= p̄(x)+ p*(x, y) (7) where p̄ varies linearly in the flow direction (i.e. p̄ =bx, where b is a constant), while p*(x, y) is a periodic part of the pressure, which repeats itself in each modulus. Eqs. (2), (3) and (5), plus the continuity equation (Eq. (1)) complete the problem formulation, together with appropriate periodically fully developed (cyclic) boundary conditions. These boundary conditions are A.O. Nieckele et al. / J. Non-Newtonian Fluid Mech. 75 (1998) 43–54 46 the symmetry condition at the longitudinal boundaries, which imply zero normal derivatives for all variables, except for the normal velocity component, which is null; “ the periodically fully developed flow condition at the inlet and outlet boundaries; “ The no-slip and impermeability conditions at the tube walls, which imply zero velocity. Assuming two-dimensional flow, the governing equations can be written in the Cartesian coordinate system shown in Fig. 1, and made dimensionless through the use of the following dimensionless variables: “ u U= , ū P( = V= p̄ , rū 2 x x*= , D 6 ū P= y*= (8) p* , rū 2 ū = m; 1 rbL (1− f) y D (9) (10) where u and 6 are, respectively, the x- and y-components of the velocity vector, m; the mass flow rate between two rows of the tube bundle, L the row spacing in the y-direction (Fig. 1), D the tube outer diameter (Fig. 1), b the tube length in the z-direction and f is the porosity, given by f= 1− p D2 8 SL (11) The physical dimensionless parameters that govern the problem are the Reynolds number, Re, and the yield number, Y, defined respectively as Re= rūD , mp Y= t0 mpū/D (12) In addition, there are two dimensionless geometrical parameters, namely, the ratios D D*= , S L*= L S (13) where S is the column spacing (Fig. 1). Some of the pressure drop results will be presented in the form of the product fRe, where f is the friction factor, defined below. f f3 f3 Dp̄ D bD = 3 3 (1 −f ) 1 2 S (1 − f ) 1 2 rū rū 2 2 (14) Clearly, f is a dimensionless form of the overall pressure drop across the computational domain, Dp̄, which is equal to bS. A.O. Nieckele et al. / J. Non-Newtonian Fluid Mech. 75 (1998) 43–54 47 Fig. 2. The mesh. 3. Numerical method The conservation equations were discretized with the aid of the finite volume method described in [18], using the central difference scheme. A non-orthogonal curvilinear system of coordinates which adapts to the boundaries of the domain was employed. Staggered velocity components were used to avoid unrealistic pressure fields, and the contra-variant velocity component was selected as the dependent variable in the momentum conservation equations [19,20]. The pressure-velocity coupling was solved by an algorithm based on SIMPLEC [21]. the resulting algebraic system was solved via the TDMA line-by-line algorithm [18] with the block correction algorithm [22] to increase the convergence rate. A transfinite interpolation scheme was employed to generate a mesh with 80×40 control volumes. A typical mesh is illustrated in Fig. 2. Extensive mesh tests were performed in order to assure essentially mesh-independent results. To this end, results were also obtained with a 128 × 60 mesh, for a Newtonian fluid and for a viscoplastic material with Y =100. Comparisons between values for the product between friction factor and Reynolds number, fRe, obtained with the 80 × 40 and 128 × 60 meshes showed a 1% difference for the Newtonian case and a 1.66% difference for the viscoplastic case. The fRe result obtained for a Newtonian fluid with D*= 0.66, L*= 0.58 and a 80×40 mesh, namely, fRe= 136.8, was also compared with the fRe value calculated from the data presented in [1], which turned out to be fRe= 150. Therefore, the agreement with experimental data was within 9% for this case. In order to check the implementation of the GNL with a biviscosity function, the flow through a straight tube was also simulated, and an excellent agreement with the Bingham plastic analytical solution was obtained. For example, for Y =10 and a mesh with 42 control volumes in the radial direction, analytical and numerical axial velocity profiles were graphically indistinguishable, and the error in fRe was equal to 0.4%. 4. Results and discussion Because most viscoplastic materials are highly viscous, in the present paper we focus our attention on flows with negligible inertial forces. To this end, although inertia terms are present 48 A.O. Nieckele et al. / J. Non-Newtonian Fluid Mech. 75 (1998) 43–54 Fig. 3. Velocity field, Y =5. Fig. 4. Velocity field, Y =500. in the formulation and in the numerical code employed, trial runs showed that when Re 51 inertial effects were always negligible. Therefore, the Reynolds number was kept constant and equal to unity for all cases. Velocity fields are shown in Figs. 3 and 4 for Y =5 and Y =500. Firstly, it can be seen that the flow field smoothly contours the rods for both values of Y. In the throat region the velocity Fig. 5. X-component velocity profile at throat. A.O. Nieckele et al. / J. Non-Newtonian Fluid Mech. 75 (1998) 43–54 49 Fig. 6. Viscosity field, Y =5. Fig. 7. Viscosity field, Y= 100. gradient near the tube wall is mild for low values of Y, and sharp for high Y’s. Away from the tube wall, where shear stresses are low, the material flows with essentially uniform velocity, due to the yield stress effect. This low-deformation-rate region is larger for the large Y case, as expected. The region of low deformation rate at the throat is further illustrated in Fig. 5. In this figure, the velocity profile of the x-component of velocity, u, is plotted as a function of y along the throat, for Y = 5 and 500. The mild variation of the velocity in the region of low deformation rate would not exist for a true Bingham plastic, but is consistent with the biviscosity law. Regions of high viscosity (and hence low deformation rate) are visualized in Figs. 6–8. In these figures, the viscosity increases with darkness, and, in the black regions, it is highest and equal to m . In the following discussion, these black regions are called high-viscosity regions. Fig. 8. Viscosity field, Y= 500. 50 A.O. Nieckele et al. / J. Non-Newtonian Fluid Mech. 75 (1998) 43–54 Fig. 9. Lines of constant periodic pressure, p*. Y= 5. Fig. 10. Lines of constant periodic pressure, p*. Y= 500. For Y = 5 (Fig. 6), it is seen that only at the throat and away from the rod a high viscosity region is observed. For higher yield stresses, however (Figs. 7 and 8), two large high-viscosity regions appear, one in the throat region and the other in the neighborhood of the domain midplane. For Y =5, the width of the high-viscosity region at the throat is about half the gap, Fig. 11. Pressure drop vs D*. A.O. Nieckele et al. / J. Non-Newtonian Fluid Mech. 75 (1998) 43–54 51 Fig. 12. Pressure drop vs L*. Fig. 13. Pressure drop vs Y. while for Y = 100 it already occupies most of the gap. The velocity and viscosity gradients are quite high a the throat and near the wall for Y= 100 and 500. Experimentally, it would probably be easy to mistakenly interpret such sharp gradients as wall slip. Moreover, for Y = 500 (Fig. 8), small high-viscosity regions are observed close to the rods. In the sharply accelerating and decelerating portions of the flow there is no high-viscosity region, not even for the largest Y case (Fig. 8). Despite the relatively low shear stresses found in these regions, the normal extra-stress components are responsible for large extra-stress magnitudes, which exceed the yield stress. It is illustrative to observe, in Fig. 4, the velocity field at the high-viscosity regions shown in Fig. 8. It becomes evident that these high-viscosity regions are not exactly rigid or plug-flow regions, because throughout them the velocity varies slightly, especially in direction. True Bingham plastics are perfectly rigid at low stresses, but small deformation rates at high-viscosity regions are in agreement with the biviscosity law. Figs. 9 and 10 illustrate the cyclic pressure distribution for Y=5 and 500, respectively. It is seen that sharp pressure variations occur near the tube wall for the large Y case, as expected. A.O. Nieckele et al. / J. Non-Newtonian Fluid Mech. 75 (1998) 43–54 52 Table 1 Newtonian values of fRe fReN D* = 0.47 L* =0.5 D* =0.65 L* = 0.7 D* =0.56 L*= 0.5 D* = 0.65 L* =0.58 D*= 0.56 L* = 0.38 189.95 132.73 161.35 136.8 361.0 Moreover, the x-dependence of the periodic pressure is such that it falls along the throat regions, and rises with x elsewhere. this is readily understood with the aid of Eq. (7), which is essentially the definition of p*. By combining the just described trend for p* with the linear term −bx, we obtain a total pressure that is essentially constant away from the rods, and that falls sharply in their neighborhood. This is a physically plausible trend, since it is expected that viscous dissipation of mechanical energy be confined to the rod vicinity. Pressure drop results are presented in Figs. 11–13 in the form bN/b, where bN is the value of b relative to the corresponding Newtonian case. Some results of Newtonian fRe values are presented in Table 1. For a given Newtonian fluid, flow rate and geometry, the corresponding bN value can be obtained if fRe is known. From Fig. 11 it can be observed that, as the tube diameter is increased (at constant flow rate), pressure drop is decreased. This trend is explained as follows. As the tube diameter is increased, the passages become narrower, and the deformation rates become larger. From Eq. (5) it is seen that, as g; becomes large, the importance of the yield stress diminishes, the viscosity decreases, and its behavior approaches the Newtonian behavior. On the other hand, increasing the row spacing, L, causes pressure drop to increase, as it can be observed from Fig. 12. This is so because larger row spacings imply wider passages, and hence lower deformation rates. Again, from Eq. (5) it is clear that lower deformation rates yield larger viscosities and hence larger pressure drops. Fig. 13 shows that increasing the yield stress causes pressure drop to increase, because increasing the yield stress gives larger viscosities. It can be seen that, for Y =500, pressure drop is two or three orders of magnitude larger than the corresponding Newtonian value, depending on geometry. Fig. 14. Flow rate vs pressure drop. A.O. Nieckele et al. / J. Non-Newtonian Fluid Mech. 75 (1998) 43–54 53 Finally, Fig. 14 gives the dimensionless flow rate, Q* =m; mp/rbt0L 2, as a function of the dimensionless pressure drop, b* = bS/t0. This figure clearly illustrates that there is a minimum pressure drop below which there is no flow. These critical values of b* correspond to the Y “ limits for the different geometries. 5. Summary and conclusions The flow of viscoplastic materials through arrays of tube bundles in staggered arrangement is studied. Velocity, viscosity and pressure fields are presented and discussed. Pressure drop results are also presented as a function of the Yield number and geometrical parameters. It is shown that, near the throats formed between adjacent tubes, a high-viscosity region appears in the central core, away from the tube walls. This low-deformation-rate region occurs due to the low extra-stress level which does not exceed the material’s yield stress. Due to the biviscosity model employed, which is slightly different from the Bingham plastic model, the high-viscosity regions observed do not have exactly uniform velocity distributions. Changing the geometrical parameters such that the passages become narrower causes pressure drop to decrease, because the deformation rates increase, decreasing the viscosity of the material. The larger is the material’s yield stress, the larger is the pressure drop, because, as the yield stress is increased, the viscosity function increases its value at a given deformation rate. There is a critical pressure drop below which there is no flow, in the limit Y “ . This critical value depends on tube spacing and arrangement. Acknowledgements This research program has been executed with financial support of Petrobras S.A., CNPq and FAPERJ. References [1] C. Chmielewski, K. Jayaraman, The effect of polymer extensibility on crossflow of polymer solutions through cylinder arrays, J. Rheol. 36 (1992) 1105– 1126. [2] A. Z& ukauskas, Heat transfer from tubes in crossflow, in: Advances in Heat Transfer, Vol. 8, Academic Press, New York, 1972, pp. 93 – 160. [3] T. Mizushina, H. Usui, Approximate solution of the boundary layer equations for the flow of a non-Newtonian fluid around a circular cylinder, Heat Trans. Jpn. Res. 7 (1978) 83 – 92. [4] T. Mizushina, H. Usui, K. Veno, T. Kato, Experiments of pseudoplastic fluid cross flow around a circular cylinder, Heat Trans. Jpn. Res. 7 (1978) 92 – 101. [5] C. Chmielewski, C.A. Petty, K. Jayaraman, Crossflow of elastic liquids through arrays of cylinders, J. Non-Newtonian Fluid Mech. 35 (1990) 309 – 325. [6] R.B. Bird, G.C. Dai, B.J. Yarusso, The rheology of flows of viscoplastic materials, Rev. Chem. Eng. 1 (1983) 1 – 70. [7] G.G. Lipscomb, M.M. Denn, Flow of Bingham fluids in complex geometries, J. Non-Newtonian. Fluid Mech. 14 (1984) 337 – 346. 54 A.O. Nieckele et al. / J. Non-Newtonian Fluid Mech. 75 (1998) 43–54 [8] D.K. Gartling, N. Phan-Thien, A numerical simulation of a plastic fluid in s parallel-plate plastometer, J. Non-Newtonian Fluid Mech. 14 (1984) 347 – 360. [9] E.J. O’Donovan, R.I. Tanner, Numerical study of the Bingham squeeze film problem, J. Non-Newtonian Fluid Mech. 15 (1984) 75 – 83. [10] T.C. Papanastasiou, Flows of materials with yield, J. Rheol. 31 (1987) 385 – 404. [11] K.R.J. Ellwood, G.C. Georgiou, T.C. Papanastasious, J.O. Wilkes, Laminar jets of Bingham-plastic liquids, J. Rheol. 34 (1990) 787 – 812. [12] S.S. Abdali, E. Mitsoulis, N.C. Markatos, Entry and exit flows of Bingham fluids, J. Rheol. 36 (1992) 389. [13] C.R. Beverly, R.I. Tanner, Numerical analysis of three-dimensional Bingham plastic flow, J. Non-Newtonian Fluid Mech. 42 (1992) 85 – 115. [14] S.D.R. Wilson, Squeezing flow of a yield-stress fluid in a wedge of slowly-varying angle, J. Non-Newtonian Fluid Mech. 50 (1993) 45 – 63. [15] S.D.R. Wilson, A.j. Taylor, the channel entry problem for a yield stress fluid, J. Non-Newtonian Fluid Mech. 65 (1996) 165 – 176. [16] J.M. Piau, Flow of a yield stress fluid in a long domain. Application to flow on an inclined plane, J. Rheol. 40 (1996) 711 – 723. [17] S.V. Patankar, C.H. Liu, E.M. Sparrow, Fully developed flow and heat transfer in ducts having streamwise periodic variations of cross-sectional area, ASME J. Heat Trans. 99 (1977) 180 – 186. [18] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, New York, 1980. [19] L.F.G. Pires, A numerical method for the solution of flows using contravariant components in non-orthogonal coordinates, PhD Thesis (in Portuguese), Pontifı́cia Universidade Católica-RJ, 1995. [20] L.F.G. Pires, A.O. Nieckele, A numerical method for the solution of flows using contravariant components in non-orthogonal coordinates, Proc. V Brazillian Meet. on Thermal Sciences, SP (in Portuguese), 1994, pp. 343 – 346. [21] J.P. Van Doormaal, G.D. Raithby, Enhancements of the SIMPLE method for prediction incompressible fluid flows, Num. Heat Trans. 7 (1984) 147 – 163. [22] A. Settari, K. Aziz, A generalization of the additive correction methods for the iterative solution of matrix equations, SIAM J. Num. Anal. 10 (1973) 506 – 521. .