- Email: [email protected]
Interactions of two rigid spheres translating collinearly in creeping flow in a Bingham material
J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
Interactions of two rigid spheres translating collinearly in creeping flow in a Bingham material Benjamin T. Liu a,b , Susan J. Muller a,b,∗ , Morton M. Denn c a
c
Materials Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720-1462, USA b Department of Chemical Engineering, University of California, Berkeley, CA 94720-1462, USA The Benjamin Levich Institute for Physico-Chemical Hydrodynamics and Department of Chemical Engineering, City College of the City University of New York, New York, NY 10031, USA Received 15 January 2003; received in revised form 25 April 2003
Abstract The interactions of two identical rigid spheres of radius R translating in a Bingham material in creeping flow along their line of centers are calculated at various sphere separations using the finite element method. The yield surfaces are determined by an extrapolation using a regularized constitutive model. Two spheres falling in a line interact at separations greater than that corresponding to the superposition of yield surfaces for isolated spheres falling in an unbounded medium. Three distinct regimes are identified: for L/R > 5.5, the spheres move in separate yield envelopes; for 4 < L/R < 5.5, the spheres fall in a single coalesced envelope in which at least part of the region between the spheres is yielded; for L/R < 4 the spheres fall in a coalesced envelope and are connected by an unyielded plug. Drag reduction of up to 30% relative to the single-sphere case is observed as the sphere separation is decreased. The yield surfaces for two approaching spheres show a shorter range of interaction compared to the case of two falling spheres; a very slight drag reduction relative to the single-sphere case is observed for two approaching spheres. © 2003 Elsevier B.V. All rights reserved. Keywords: Yield stress; Bingham material; Regularization; Viscoplasticity; Particle interactions
1. Introduction Flows of yield stress, or viscoplastic, materials can be important in a variety of applications. Such materials deform as non-Newtonian viscous liquids when the local stress exceeds a critical value; at lower levels of stress they behave as solids. Few two-dimensional flows of these materials have been addressed in the literature because of the complexity introduced by the discontinuous constitutive description. One notable exception is the problem of creeping flow past a rigid sphere, either in an unbounded regime [1–3] or along the axis of a cylinder [4,5]. In this work, we look at the interactions of two rigid spheres ∗
Corresponding author. Tel.: +1-510-642-4525; fax: +1-510-642-4778. E-mail address: [email protected] (S.J. Muller). 0377-0257/$ – see front matter © 2003 Elsevier B.V. All rights reserved. doi:10.1016/S0377-0257(03)00111-3
50
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
moving along their line of centers. Such interactions can be important both in the settling of solid particles and the transport, coalescence, and release of gas bubbles. Additionally, by incrementally increasing the complexity of the system, we seek to identify computational problems associated with viscoplastic flows that are not apparent in simpler systems. A model describing the deviatoric stress tensor (τ ) in a yield stress material, employing a von Mises yield criterion, was formulated by Oldroyd [6,7] and Prager [8], as follows: 1 τY τ = ηp + √ D for IIτ > τY , (1) 2 1/2(IID ) 1 II ≤ τY , (2) D = 0 for 2 τ where τ Y and ηp are the yield stress and the “plastic viscosity”, respectively; D is the rate of deformation tensor, defined in terms of the velocity v as D = ∇v + (∇v)T ; the second invariants IID and IIτ are defined D:D and τ :τ , respectively. This material is known as a Bingham material. Flows of Bingham materials can be characterized by a Bingham number, defined 2τY R , (3) Bn = ηp V where R and √ V are characteristic length and velocity scales, respectively. Regions that satisfy the yield criterion ( 1/2(IIτ ) > τY ) are referred to as yielded, while those that do not satisfy the yield criterion are unyielded and assumed to behave as rigid solids. The boundaries between yielded and unyielded regions are known as yield surfaces. The discontinuity in the Bingham equation creates significant difficulties in the numerical solution of viscoplastic flow problems. To address this problem, many authors have used regularized constitutive equations, in which the discontinuous Bingham equation is replaced by a continuous equation characterized by a regularization parameter; for two-dimensional creeping flow, the flow fields calculated from these models converge to the Bingham solution as the regularization parameter tends to infinity [9]. A regularized constitutive equation is used in the present study and is defined in Section 3. 2. Previous work Several authors have addressed the problem of a single sphere translating in a Bingham material. Most of these studies have focused on determining the location of the yield surfaces by using regularized constitutive equations. Beris et al. [2] formulated the problem as a free-boundary problem, where the yield surface was treated as an unknown boundary and the calculation was restricted to the yielded region. While this work is considered to be the benchmark for the single-sphere problem, the approach requires a priori knowledge of the general location and number of yield surfaces. A more recent study by Liu et al. [1] identified the positions and shapes of the yield surfaces by performing the calculation over a fixed domain. The flow field was simulated in both the yielded and unyielded regions, so no prior knowledge of the yield surfaces was required. It was found that the shape and extent of the yielded region was dependent on the value of the regularization parameter; moreover,
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
51
the limit of very large regularization parameter could not be reached without numerical error. A method of extrapolating the location of the yield surfaces using a sequence of values of the regularization parameter was developed. This method, which was used in the present study, is described in Section 3. Beris et al., Liu et al. and other researchers have used variational integrals to evaluate the accuracy of their numerical solutions. Prager [10] defines such variational integrals for creeping flow of a Bingham material in a fluid volume Ω, bounded by a surface SΩ with a normal n on which the velocity is fixed, as follows: 1 1 1 (4) ηp IID + 2τY IID dΩ, H= πηp RV2 Ω 2 2
2
1 1 1 1
− K= IIτ − τY + IIτ − τY dΩ + 2 n · [τ · v] dSΩ .
πηp RV2 4 Ω 2 2 SΩ
(5)
For the exact solution, H and K take on the same numerical value. These variational integrals are also used in the present work to verify the accuracy of numerical solutions. The interactions of pairs of particles in Newtonian flows have been extensively studied. Creeping flow of two spheres falling along the line of their centers in an infinite Newtonian fluid has been solved exactly by Stimson and Jeffery [11]. An analytical solution has also been obtained for two spheres approaching one another in creeping flow of a Newtonian fluid [12]. The spheres in a creeping Newtonian fluid exhibit long-range interactions; these interactions cause a decrease in the drag on each sphere for spheres falling in a line, whereas the drag on each sphere increases for approaching spheres.
3. Formulation We consider the problem of two identical rigid spheres of radius R separated by a distance L translating along the axis connecting their centers in an incompressible, isothermal Bingham material of infinite extent. We assume creeping flow. The flow is assumed to be axisymmetric about the axis of translation, so we use a cylindrical (r, z) coordinate system centered on one of the spheres, as shown in Fig. 1. We impose the axisymmetry condition (vr = 0, Fz = 0) along the central axis (r = 0). Both velocity components vanish in the far-field (vr = vz = 0). For the case of two spheres falling in a line, the velocity on the surface of each sphere is set to V in the negative z-direction (vr = 0, vz = −V ). From the creeping flow assumption, we have fore-aft symmetry around the two-sphere system, so we can enforce a fore-aft plane of symmetry (vr = 0, Fz = 0, where Fz is the z-component of the boundary traction) along the mid-plane between the spheres (z = −L/2). The fore-aft symmetry (and hence reversibility) of the flow can be assumed despite the non-linearity of the constitutive relation because the apparent viscosity is an even function of the velocity. This reversibility also implies that the spheres will not exert a net attractive/repulsive force on each other, so both spheres fall at the same fixed velocity V and the flow is steady. For the case of two spheres approaching one another, we set the velocities on the spheres to be equal and opposite. The creeping flow assumption results in mirror symmetry about the mid-plane between the spheres, so we set the velocity on the surface of the first sphere as before (vr = 0, vz = −V ), but enforce a plane of mirror symmetry (vz = 0, Fr = 0, where Fr is the r-component of the boundary traction) along
52
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
z
z
r
R
r
R
V
V L
L fore-aft symmetry plane
mirror symmetry plane
V R
R V
(a)
(b)
Fig. 1. Schematic representation of (a) two spheres falling collinearly and (b) two approaching spheres.
the mid-plane between the spheres (z = −L/2). In this case, the sphere separation is decreasing and the flow is inherently unsteady, but the time dependent terms can be neglected in creeping flow. Various values of sphere separation, L, were considered. We solved the equations of mass and momentum conservation using a penalty finite element method [13]. We used a regularized constitutive equation [1,2,14], in which the dimensionless stress and velocity are related as follows:
(1/2)Bn τ = 1+ √ D. (6) (1/2)IID + PB−1 PB is a dimensionless regularization parameter and Bn is defined using R as the characteristic length and V as the characteristic velocity. The details of the calculation method are outlined in our previous study [1]. The following convergence criteria was used for the non-linear iterations: (v i − v i−1 ) · (v i − v i−1 ) < 10−6 , (7) vi · vi where v i is the vector of nodal velocities from the current iteration and v i−1 is the vector of nodal velocities from the previous iteration; this criterion was also used in our earlier work but reported incorrectly [1]. The value of the penalty parameter used here (107 ) is also identical to that used in our previous study. We use Bn = 340.7 unless otherwise noted. This relatively high value was chosen in order to focus on
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
53
Fig. 2. Representative mesh used for calculations. Mesh F1 from Table 1.
the yield stress contributions and corresponds to a value used in previous studies [1,2]. We approximate the infinite medium with a finite mesh domain characterized by Rmesh and Zmesh . Fig. 2 shows the domain and refinement for a characteristic mesh; the meshes used for the simulations are summarized in Table 1. The mesh refinement was chosen such that the velocities for Bn = 340.7 were indistinguishable when the number of elements was doubled; the variational integral H differed by less than 0.2% with mesh doubling, while K was within 1.5%. The slight errors in the strain rate field due to inadequate mesh refinement are visible in the strain rate contours (shown in Fig. 3); the finer meshes were used to ensure that the strain rates were smooth. The velocity and strain rate fields for several different values of PB were used in the extrapolation to define the yield surfaces. To facilitate this calculation, we define a scaled strain rate by IID ≡
IID , IID,Y
(8)
54
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
Table 1 Refinement and characteristic lengths of meshes used for the simulations Mesh
L/R
Rmesh
Zmesh
Elements
Nodes
A1 A2 B C D E F1 F2 G H
7.0 7.0 6.0 5.5 5.0 4.5 4.0 4.0 2.5 2.1
5.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0
3.5 10.0 10.0 10.0 10.0 10.0 10.0 5.0 10.0 10.0
10100 5564 27434 24990 27397 24043 20725 9600 11698 10238
40871 22537 110313 125531 110165 96729 83435 38881 58939 51639
Meshes A2 and F2 have been selectively refined near the spheres relative to meshes A1 and F1, respectively.
where IID,Y is the strain rate at yield, obtained by solving Eq. (6) at the yield stress: −1 2 1 IID,Y = 2 (−PB + PB−2 + 2BnP−1 B ) .
(9)
The scaled strain rate can be thought of as the relative state of yield; IID = 1 corresponds to the yield surfaces defined by the calculated stress for a fixed value of PB . As we increase the value of PB , the position of the yield surface changes; the true yield surfaces are found in the Bingham limit (PB → ∞). We note that IID increases in some regions while it decreases in others as we increase PB ; a representative profile is shown in Fig. 4. Moreover, this behavior is observed up to values of PB beyond which numerical results are unreliable. Extrapolating to the Bingham limit, we identify the regions where IID increases with PB as yielded, and the regions where IID tends to zero with increasing PB as unyielded. We can
Fig. 3. Contours of strain rate for two spheres falling collinearly calculated using (a) a coarse mesh (F2) and (b) a fine mesh (F1). L/R = 4.0, Bn = 340.7, PB = 5000. The contours range logarithmically from 10−14 to 106 , in 20 intervals.
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
55
Fig. 4. Profiles of the scaled strain rate (IID ) vs. radial position along the mid-plane (z = −L/2) for various values of PB . Two spheres falling collinearly: L/R = 4.0 and Bn = 340.7.
56
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
then define the extrapolated yield surfaces as the loci of points where IID is invariant with respect to PB . We can use this extrapolation to identify yielded regions even where the calculated scaled strain rates are less than unity (as in Fig. 4b); the sensitivity of the results to PB must be monitored, however, to ensure that the regions of increasing and decreasing IID are indeed tending to infinity and zero, respectively. For each value of L/R, simulations were run for different values of PB . We define a quantity φ≡
∂ log IID , ∂ log PB
(10)
which represents the change in relative state of yield with regularization parameter. The value of φ was calculated using a linear least squares fit at each nodal point and interpolated within the mesh. The extrapolated yield surfaces were defined by φ = 0; regions of positive and negative φ were identified as yielded and unyielded, respectively. Preliminary studies were performed using PB = 5 × 103 , 6 × 103 , 7 × 103 , 8 × 103 , 9 × 103 , and 104 ; the least squares fit using the linear approximation was very good (the coefficient of determination, r2 , was found to be greater than 0.95 over nearly the entire field). Moreover, the results using the two values PB = 5 × 103 and 8 × 103 were identical to the results using the full range of PB (5 × 103 , 6 × 103 , 7 × 103 , 8 × 103 , 9 × 103 , 104 ). We thus used PB = {5 × 103 , 8 × 103 } for the present calculations. Simulations using larger values of PB (up to 106 ) also showed good agreement, though numerical round-off error was readily apparent in some of these cases. The Stokes drag coefficient, defined CS =
Fdrag , 6πηp VR
(11)
where Fdrag is the drag experienced by each sphere, was also calculated for varying sphere separations. 4. Results: two spheres falling collinearly 4.1. Yield surfaces In calculating the locus of points where φ = 0, several different surfaces were identified, some of which may represent true extrapolated yield surfaces and some of which are mesh related artifacts. Fig. 5a shows the various surfaces for two spheres falling in a line with L/R = 7.0. As expected, we see an envelope surrounding each sphere (α), corresponding to the yield envelope observed in previous studies [2,4]. The yielded envelope is nearly identical to the envelope calculated for a single sphere in our previous work [1]. In addition, we see apparent yield surfaces near the corners of the mesh (β, γ) and also near the mid-plane between the spheres (δ). When we plot the values of the strain rate along the apparent yield surfaces, as in Fig. 5b, we find that the strain rates associated with surfaces β, γ, and δ, are several orders of magnitude lower than those on the yield envelope. Moreover, regions β and γ were dependent on the size of the computational domain (this behavior was also observed for single-sphere calculations), and the regions moved when the calculation was repeated on a larger mesh (Rmesh = Zmesh = 10). Region δ disappeared with selective mesh refinement. Examination of the strain rates near the mid-plane (Fig. 6) show that
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
57
Fig. 5. (a) Calculated yield surface for two spheres falling collinearly: Bn = 340.7 and L/R = 7.0, using mesh A1. Unyielded regions are shaded. β, γ, δ indicate anomalous yield surfaces. (b) Strain rate plotted against radial position along calculated yield surfaces; two spheres falling collinearly: Bn = 340.7 and L/R = 7.0.
58
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
Fig. 6. Profiles of the scaled strain rate (IID ) vs. radial position along the mid-plane (z = −L/2). Coarse mesh (A1) results in spurious yield surfaces. Two spheres falling collinearly: Bn = 340.7 and L/R = 7.0.
this apparent yielded region (δ) actually corresponds to the strain rate passing through zero; because this results in a cusp in IID , small mesh-dependent errors can cause the sign of φ to change, resulting in an apparent yield surface when the mesh is insufficiently refined. Hence, we clearly see that surfaces β, γ, and δ of Fig. 5a are mesh-related artifacts, and increased local mesh refinement can be used to eliminate these regions (the remaining meshes were refined to eliminate or minimize the mid-plane regions). We can also use the magnitude of the strain rate along the calculated yield surface to discriminate the true yield surfaces from these spurious mesh-related yield surfaces. This method is particularly useful for discriminating the corner regions, as it is computationally inefficient to refine the mesh regions far away from the spheres. Fig. 7 shows the true yield surfaces for two spheres falling in a line at smaller values of L/R (only the portion of the computational region near the yielded region is shown in this and subsequent figures). These surfaces were determined from the calculated yield surfaces by eliminating those surfaces where IID is less than some threshold value, as determined by a plot analogous to Fig. 5b. This threshold value is dependent on the local refinement of the mesh, and was found to be around 10−7 for all of the meshes except meshes G and H, where 10−5 was used. At L/R = 6.0 we see that the yield
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
59
Fig. 7. Calculated yield surfaces for two spheres falling collinearly: Bn = 340.7 and (a) L/R = 6.0; (b) L/R = 5.0; (c) L/R = 4.5; (d) L/R = 4.0; (e) L/R = 2.5; (f) L/R = 2.1. Spurious mesh-related yield surfaces have been removed.
envelopes are separate and essentially unchanged from those calculated for a single sphere. For these larger separations (L/R ≥ 6.0), there is essentially no interaction between the spheres. We can see that as we decrease the sphere separation, the spheres begin to interact and the yield envelopes eventually coalesce.
60
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
Fig. 8. Streamlines of the re-circulation pattern within the yielded envelope, two spheres falling collinearly: Bn = 340.7 and L/R = 4.0; 20 evenly spaced streamlines (ψ = −0.03–0.5). The axial velocity (vz ) profile along the mid-plane (z = −2.0) is shown for reference.
For L/R ≤ 5.0, the two spheres are no longer contained in separate yield envelopes, but rather in a single coalesced envelope. Although the yield surfaces for two isolated spheres would not intersect at L/R = 5, the stress field between the two spheres is not a linear superposition of two isolated spheres, so it is not surprising that the spheres are found to interact at this distance. We note distinct unyielded regions within the yielded envelope: along the mid-plane, centered at r ∼ 2; along the axis of symmetry, centered at the mid-plane between the spheres; and on the poles of the spheres. The region along the mid-plane corresponds to a plug-like flow of the material due to re-circulation within the yielded envelope. The re-circulation can be seen by examining the streamlines in Fig. 8 and is clearly required to preserve mass conservation within the yielded envelope. Note that the unyielded region is not in rotation, which would be inadmissible given the axisymmetric geometry, but rather moves uniformly opposite the motion of the spheres. This region must also exist as at least a point within the separate yield envelopes (L/R ≥ 6.0), although it may not be large enough to be visible in these calculations. The mid-plane unyielded region decreases in size with decreasing sphere separation and vanishes as the separation approaches 2R.
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
61
Fig. 9. Comparison of calculated yield surfaces for two spheres falling collinearly (solid lines) and a capped-cylinder (dashed lines): Bn = 340.7 and (a) L/R = 2.1; (b) L/R = 2.5; (c) L/R = 4.0; (d) L/R = 5.0.
The size of the coalesced envelope decreases with decreasing sphere separation while maintaining the same approximate shape and aspect ratio. By comparing the results with results for two spheres connected by a rigid bridge (i.e. a capped cylinder), as shown in Fig. 9, we see that the shape of the outer yield envelope is relatively insensitive to the shapes of the interior unyielded regions and appears to be dominated by the shape of the leading (and trailing) edge to the flow.
62
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
Fig. 10. Calculated yield surface for two spheres falling collinearly: Bn = 340.7 and L/R = 5.5. The annular region is expected to vanish in the Bingham limit.
The unyielded region along the axis of symmetry decreases in size with decreasing sphere separation, maintaining an approximately constant distance from the spheres, until the separation reaches about 4R; the region then becomes larger and connects the two spheres. The plug of unyielded material connecting the spheres then travels at the same velocity as the spheres. The polar caps probably exist on both the fore and aft surfaces of each sphere because of the stagnation points [2], but they are too small to be distinguishable in some of the simulations. The region on the pole closest to the mid-plane appears to increase in size with decreasing sphere separation until it abruptly coalesces when the spheres become connected at L/R = 4.0. The region on the opposite pole also appears to be increasing in size with decreasing sphere separation, but the regions are so small that the current calculations cannot resolve them accurately. The sphere separation of 5.5 (Fig. 10) represents a transition between separate yielded envelopes and a single, coalesced yielded envelope, in that the calculated yield surface appears as two separate envelopes connected by a thin annular channel. The axial velocities in this annular region and the region within it are in the same direction, however, and this configuration cannot conserve mass except through movement in the opposite direction in the surrounding unyielded region. Since the far-field velocities are fixed at zero, such a flow cannot occur in the Bingham limit, and the calculated velocities in the apparent plug do, in fact, decrease with increasing PB , as does the thickness of the annulus. A regularization parameter sufficiently large to see the complete vanishing of the annular channel was inaccessible for this transition case because of round-off error. The yield surfaces for smaller values of Bn were calculated for L/R = 4.0 and are shown in Fig. 11. The yield surfaces are relatively insensitive to Bn in this range; the height of the yielded envelope along the axis of symmetry (measured from the center of the nearest sphere) increases by only 14%, and the width at the mid-plane increases by only 11%, as Bn is decreased from 340.7 to 14.91. These results are similar to previous single-sphere calculations, which show slightly larger increases (20–30%) over the same range of Bn [1,2].
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
63
Fig. 11. Calculated yield surfaces for two spheres falling collinearly with L/R = 4.0 at various values of Bn.
4.2. Drag coefficients and variational integrals The calculated values of CS , H and K are insensitive to the value of PB ; they vary by less than 0.5% with a twofold change in the regularization parameter (5 × 103 –104 ). The accuracy of the results is confirmed by the small relative difference between H and K shown in Table 2. This is true even for the calculations at L/R = 5.5, where the anomalous annular channel connecting the yield surfaces appears (cf. Fig. 10), Table 2 Calculated variational integrals for two identical rigid spheres falling along the line of centers: Bn = 340.7 and PB = 104 L/R
Mesh
H
K
(H − K)/K (%)
7.0 7.0 6.0 5.5 5.0 4.5 4.0 4.0 2.5 2.1
A1 A2 B C D E F1 F2 G H
5020.7 5075.0 5032.3 5044.5 5014.5 4915.1 4729.4 4737.7 4037.5 3820.7
5016.9 4987.2 5000.1 5003.0 4975.7 4866.5 4635.9 4584.9 3864.1 3652.2
0.08 1.8 0.64 0.83 0.78 1.0 2.0 3.3 4.5 4.6
64
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67 1 0.9 Newtonian result
normalized drag coefficient (C S /CS,single sphere)
0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 2
3
4
5
6
7
8
9
10
sphere separation (L/R)
Fig. 12. Plot of normalized drag coefficient (CS /CS,single sphere ) for two spheres falling collinearly as a function of sphere separation, Bn = 340.7. Open symbols were calculated using less refined meshes.
and underscores the sensitivity of the location of the yield surfaces to the details of the strain field. The drag coefficients normalized by the drag coefficient for a single, isolated sphere are shown for various sphere separations in Fig. 12, together with the results for a Newtonian fluid [11,15]. For L/R ≥ 6, the drag experienced by each sphere is identical to that for an isolated sphere, as expected from the shapes of the yield surfaces (Figs. 7 and 8). The range of interaction in the yield stress material is much smaller than in the Newtonian case; the presence of the second sphere reduces the drag at small separations, but to a lesser extent than in a Newtonian fluid. An approximately 30% drag reduction is observed in the limiting case when the spheres are touching (L/R = 2.0). 5. Results: two approaching spheres Several sphere separations were examined for the case of two rigid spheres approaching at equal speed. The drag experienced by each sphere is tabulated in Table 3. There is no effect on the drag force for sphere separations greater than or equal to 4.5; this contrasts with the case of two spheres falling in the same direction, where interaction is observed up to L/R = 5. For smaller sphere separations, the spheres are increasingly affected by one another. Unlike the Newtonian fluid case, however, the spheres experience a slight drag reduction rather then an augmentation [12]; the reason for the net decrease in drag is that the interactions of the two spheres causes a decrease in the local viscosity because of the shear-thinning characteristic of the Bingham material.
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
65
Table 3 Calculated variational integrals and drag correction factor for two identical rigid spheres approaching one another at the same speed: Bn = 340.7 and PB = 5 × 103 L/R
Mesh
H
K
CS /CS,single sphere
4.5 4.0 2.5
E F1 G
5032.5 5016.1 4792.4
4999.7 4979.8 4736.6
1.00 0.99 0.94
Fig. 13. Calculated yield surfaces for two approaching spheres: Bn = 340.7 and (a) L/R = 4.5; (b) L/R = 4.0.
Fig. 14. Streamlines of the re-circulation pattern within the yielded envelope, two approaching spheres: Bn = 340.7 and L/R = 4.0; 20 evenly spaced streamlines (ψ = 0.01–0.5).
66
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
The extrapolated yield surfaces for two approaching spheres are shown in Fig. 13. At L/R = 4.5, the spheres are enclosed in separate yielded envelopes, as for a single, isolated sphere; this result is consistent with the lack of interaction observed in the drag calculations. At L/R = 4.0, the yielded envelopes are coalesced. The shape of the coalesced envelope is qualitatively different from that observed for two spheres falling in a line; these differences are highlighted by comparing the streamlines for each case (Fig. 14 for two approaching spheres and Fig. 8 for two spheres falling collinearly).
6. Conclusions As in our earlier study of a single sphere in a Bingham material [1], numerical error for calculations of two spheres in creeping flow along their line of centers precludes the use of regularization parameters sufficiently large to determine the location of the yield surfaces directly from the yield condition. The yield surfaces are therefore determined from the loci of values of the normalized second invariant of the deformation rate that are invariant with respect to the regularization parameter, together with refinement of the mesh and consideration of the magnitude of the deformation rate on the apparent yield surfaces. The minimum value of the regularization parameter needed to reflect ideal Bingham behavior is dependent on the details of the flow field. For sphere separations L/R ≥ 6, the yield surfaces are identical to those for single isolated spheres. The spheres interact at separations larger than those that would be predicted from single-sphere calculations; the yielded envelopes surrounding each sphere coalesce, and internal plug-like regions appear along the symmetry axis and along the mid-plane. A decrease in the drag coefficient of up to 30% accompanies the coalescence of the yielded envelopes. The size of the coalesced yield surface in insensitive to Bingham number. The range of interaction is shorter for two approaching spheres than for two spheres falling in a line, and the shapes of the coalesced yield surfaces are qualitatively different. Only a slight drag decrease compared to the single-sphere case is predicted, but this contrasts sharply with the drag increase that is seen in a Newtonian fluid.
Acknowledgements This research was supported in part by the Department of Energy, award number DE-AC03-76SF00098, and through computing resources provided by the National Partnership for Advanced Computational Infrastructure at the University of Michigan Center for Parallel Computing. Financial support through a National Science Foundation Graduate Research Fellowship to B.T.L. is also gratefully acknowledged.
References [1] B.T. Liu, S.J. Muller, M.M. Denn, Convergence of a regularization method for creeping flow of a Bingham material about a rigid sphere, J. Non-Newtonian Fluid Mech. 102 (2002) 179. [2] A.N. Beris, J.A. Tsamopoulos, R.C. Armstrong, R.A. Brown, Creeping motion of a sphere through a Bingham plastic, J. Fluid Mech. 158 (1985) 219. [3] N. Yoshioka, K. Adachi, H. Ishimura, On creeping flow of a viscoplastic fluid past a sphere, Kagaku Kogaku 35 (1971) 1144 (in Japanese).
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
67
[4] J. Blackery, E. Mitsoulis, Creeping motion of a sphere in tubes filled with a Bingham plastic material, J. Non-Newtonian Fluid Mech. 70 (1997) 59. [5] E. Mitsoulis, Effect of rheological properties on the drag coefficient for creeping motion around a sphere falling in a tightly-fitting tube, J. Non-Newtonian Fluid Mech. 74 (1998) 263. [6] J. Oldroyd, A rational formulation of the equations of plastic flow for a Bingham solid, Proc. Camb. Philos. Soc. 43 (1947) 100. [7] J. Oldroyd, Two-dimensional plastic flow of a Bingham solid, Proc. Camb. Philos. Soc. 43 (1947) 383. [8] W. Prager, Introduction to Continuum Mechanics, Ginn, Boston, 1961. [9] R. Glowinski, J.L. Lions, R. Tremolieres, Numerical Analysis of Variational Inequalities, North-Holland, Amsterdam, 1981. [10] W. Prager, in: G. Birkhoff, G. Kuerti, G. Szegö (Eds.), Studies in Mathematics and Mechanics, Academic Press, New York, 1954, p. 208. [11] M. Stimson, G.B. Jeffery, The motion of two spheres in a viscous fluid, Proc. R. Soc. London A111 (1926) 110. [12] H. Brenner, The slow motion of a sphere through a viscous fluid towards a plane surface, Chem. Eng. Sci. 16 (1961) 242. [13] M. Bercovier, M. Engelman, A finite-element method for the numerical solution of viscous incompressible flows, J. Comp. Phys. 30 (1979) 181. [14] M. Bercovier, M. Engelman, A finite-element method for incompressible non-Newtonian flows, J. Comp. Phys. 36 (1980) 313. [15] J. Happel, R. Pfeffer, The motion of two spheres following each other in a viscous fluid, Am. Inst. Chem. Eng. J. 6 (1960) 129.
Interactions of two rigid spheres translating collinearly in creeping flow in a Bingham material Benjamin T. Liu a,b , Susan J. Muller a,b,∗ , Morton M. Denn c a
c
Materials Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720-1462, USA b Department of Chemical Engineering, University of California, Berkeley, CA 94720-1462, USA The Benjamin Levich Institute for Physico-Chemical Hydrodynamics and Department of Chemical Engineering, City College of the City University of New York, New York, NY 10031, USA Received 15 January 2003; received in revised form 25 April 2003
Abstract The interactions of two identical rigid spheres of radius R translating in a Bingham material in creeping flow along their line of centers are calculated at various sphere separations using the finite element method. The yield surfaces are determined by an extrapolation using a regularized constitutive model. Two spheres falling in a line interact at separations greater than that corresponding to the superposition of yield surfaces for isolated spheres falling in an unbounded medium. Three distinct regimes are identified: for L/R > 5.5, the spheres move in separate yield envelopes; for 4 < L/R < 5.5, the spheres fall in a single coalesced envelope in which at least part of the region between the spheres is yielded; for L/R < 4 the spheres fall in a coalesced envelope and are connected by an unyielded plug. Drag reduction of up to 30% relative to the single-sphere case is observed as the sphere separation is decreased. The yield surfaces for two approaching spheres show a shorter range of interaction compared to the case of two falling spheres; a very slight drag reduction relative to the single-sphere case is observed for two approaching spheres. © 2003 Elsevier B.V. All rights reserved. Keywords: Yield stress; Bingham material; Regularization; Viscoplasticity; Particle interactions
1. Introduction Flows of yield stress, or viscoplastic, materials can be important in a variety of applications. Such materials deform as non-Newtonian viscous liquids when the local stress exceeds a critical value; at lower levels of stress they behave as solids. Few two-dimensional flows of these materials have been addressed in the literature because of the complexity introduced by the discontinuous constitutive description. One notable exception is the problem of creeping flow past a rigid sphere, either in an unbounded regime [1–3] or along the axis of a cylinder [4,5]. In this work, we look at the interactions of two rigid spheres ∗
Corresponding author. Tel.: +1-510-642-4525; fax: +1-510-642-4778. E-mail address: [email protected] (S.J. Muller). 0377-0257/$ – see front matter © 2003 Elsevier B.V. All rights reserved. doi:10.1016/S0377-0257(03)00111-3
50
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
moving along their line of centers. Such interactions can be important both in the settling of solid particles and the transport, coalescence, and release of gas bubbles. Additionally, by incrementally increasing the complexity of the system, we seek to identify computational problems associated with viscoplastic flows that are not apparent in simpler systems. A model describing the deviatoric stress tensor (τ ) in a yield stress material, employing a von Mises yield criterion, was formulated by Oldroyd [6,7] and Prager [8], as follows: 1 τY τ = ηp + √ D for IIτ > τY , (1) 2 1/2(IID ) 1 II ≤ τY , (2) D = 0 for 2 τ where τ Y and ηp are the yield stress and the “plastic viscosity”, respectively; D is the rate of deformation tensor, defined in terms of the velocity v as D = ∇v + (∇v)T ; the second invariants IID and IIτ are defined D:D and τ :τ , respectively. This material is known as a Bingham material. Flows of Bingham materials can be characterized by a Bingham number, defined 2τY R , (3) Bn = ηp V where R and √ V are characteristic length and velocity scales, respectively. Regions that satisfy the yield criterion ( 1/2(IIτ ) > τY ) are referred to as yielded, while those that do not satisfy the yield criterion are unyielded and assumed to behave as rigid solids. The boundaries between yielded and unyielded regions are known as yield surfaces. The discontinuity in the Bingham equation creates significant difficulties in the numerical solution of viscoplastic flow problems. To address this problem, many authors have used regularized constitutive equations, in which the discontinuous Bingham equation is replaced by a continuous equation characterized by a regularization parameter; for two-dimensional creeping flow, the flow fields calculated from these models converge to the Bingham solution as the regularization parameter tends to infinity [9]. A regularized constitutive equation is used in the present study and is defined in Section 3. 2. Previous work Several authors have addressed the problem of a single sphere translating in a Bingham material. Most of these studies have focused on determining the location of the yield surfaces by using regularized constitutive equations. Beris et al. [2] formulated the problem as a free-boundary problem, where the yield surface was treated as an unknown boundary and the calculation was restricted to the yielded region. While this work is considered to be the benchmark for the single-sphere problem, the approach requires a priori knowledge of the general location and number of yield surfaces. A more recent study by Liu et al. [1] identified the positions and shapes of the yield surfaces by performing the calculation over a fixed domain. The flow field was simulated in both the yielded and unyielded regions, so no prior knowledge of the yield surfaces was required. It was found that the shape and extent of the yielded region was dependent on the value of the regularization parameter; moreover,
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
51
the limit of very large regularization parameter could not be reached without numerical error. A method of extrapolating the location of the yield surfaces using a sequence of values of the regularization parameter was developed. This method, which was used in the present study, is described in Section 3. Beris et al., Liu et al. and other researchers have used variational integrals to evaluate the accuracy of their numerical solutions. Prager [10] defines such variational integrals for creeping flow of a Bingham material in a fluid volume Ω, bounded by a surface SΩ with a normal n on which the velocity is fixed, as follows: 1 1 1 (4) ηp IID + 2τY IID dΩ, H= πηp RV2 Ω 2 2
2
1 1 1 1
− K= IIτ − τY + IIτ − τY dΩ + 2 n · [τ · v] dSΩ .
πηp RV2 4 Ω 2 2 SΩ
(5)
For the exact solution, H and K take on the same numerical value. These variational integrals are also used in the present work to verify the accuracy of numerical solutions. The interactions of pairs of particles in Newtonian flows have been extensively studied. Creeping flow of two spheres falling along the line of their centers in an infinite Newtonian fluid has been solved exactly by Stimson and Jeffery [11]. An analytical solution has also been obtained for two spheres approaching one another in creeping flow of a Newtonian fluid [12]. The spheres in a creeping Newtonian fluid exhibit long-range interactions; these interactions cause a decrease in the drag on each sphere for spheres falling in a line, whereas the drag on each sphere increases for approaching spheres.
3. Formulation We consider the problem of two identical rigid spheres of radius R separated by a distance L translating along the axis connecting their centers in an incompressible, isothermal Bingham material of infinite extent. We assume creeping flow. The flow is assumed to be axisymmetric about the axis of translation, so we use a cylindrical (r, z) coordinate system centered on one of the spheres, as shown in Fig. 1. We impose the axisymmetry condition (vr = 0, Fz = 0) along the central axis (r = 0). Both velocity components vanish in the far-field (vr = vz = 0). For the case of two spheres falling in a line, the velocity on the surface of each sphere is set to V in the negative z-direction (vr = 0, vz = −V ). From the creeping flow assumption, we have fore-aft symmetry around the two-sphere system, so we can enforce a fore-aft plane of symmetry (vr = 0, Fz = 0, where Fz is the z-component of the boundary traction) along the mid-plane between the spheres (z = −L/2). The fore-aft symmetry (and hence reversibility) of the flow can be assumed despite the non-linearity of the constitutive relation because the apparent viscosity is an even function of the velocity. This reversibility also implies that the spheres will not exert a net attractive/repulsive force on each other, so both spheres fall at the same fixed velocity V and the flow is steady. For the case of two spheres approaching one another, we set the velocities on the spheres to be equal and opposite. The creeping flow assumption results in mirror symmetry about the mid-plane between the spheres, so we set the velocity on the surface of the first sphere as before (vr = 0, vz = −V ), but enforce a plane of mirror symmetry (vz = 0, Fr = 0, where Fr is the r-component of the boundary traction) along
52
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
z
z
r
R
r
R
V
V L
L fore-aft symmetry plane
mirror symmetry plane
V R
R V
(a)
(b)
Fig. 1. Schematic representation of (a) two spheres falling collinearly and (b) two approaching spheres.
the mid-plane between the spheres (z = −L/2). In this case, the sphere separation is decreasing and the flow is inherently unsteady, but the time dependent terms can be neglected in creeping flow. Various values of sphere separation, L, were considered. We solved the equations of mass and momentum conservation using a penalty finite element method [13]. We used a regularized constitutive equation [1,2,14], in which the dimensionless stress and velocity are related as follows:
(1/2)Bn τ = 1+ √ D. (6) (1/2)IID + PB−1 PB is a dimensionless regularization parameter and Bn is defined using R as the characteristic length and V as the characteristic velocity. The details of the calculation method are outlined in our previous study [1]. The following convergence criteria was used for the non-linear iterations: (v i − v i−1 ) · (v i − v i−1 ) < 10−6 , (7) vi · vi where v i is the vector of nodal velocities from the current iteration and v i−1 is the vector of nodal velocities from the previous iteration; this criterion was also used in our earlier work but reported incorrectly [1]. The value of the penalty parameter used here (107 ) is also identical to that used in our previous study. We use Bn = 340.7 unless otherwise noted. This relatively high value was chosen in order to focus on
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
53
Fig. 2. Representative mesh used for calculations. Mesh F1 from Table 1.
the yield stress contributions and corresponds to a value used in previous studies [1,2]. We approximate the infinite medium with a finite mesh domain characterized by Rmesh and Zmesh . Fig. 2 shows the domain and refinement for a characteristic mesh; the meshes used for the simulations are summarized in Table 1. The mesh refinement was chosen such that the velocities for Bn = 340.7 were indistinguishable when the number of elements was doubled; the variational integral H differed by less than 0.2% with mesh doubling, while K was within 1.5%. The slight errors in the strain rate field due to inadequate mesh refinement are visible in the strain rate contours (shown in Fig. 3); the finer meshes were used to ensure that the strain rates were smooth. The velocity and strain rate fields for several different values of PB were used in the extrapolation to define the yield surfaces. To facilitate this calculation, we define a scaled strain rate by IID ≡
IID , IID,Y
(8)
54
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
Table 1 Refinement and characteristic lengths of meshes used for the simulations Mesh
L/R
Rmesh
Zmesh
Elements
Nodes
A1 A2 B C D E F1 F2 G H
7.0 7.0 6.0 5.5 5.0 4.5 4.0 4.0 2.5 2.1
5.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0
3.5 10.0 10.0 10.0 10.0 10.0 10.0 5.0 10.0 10.0
10100 5564 27434 24990 27397 24043 20725 9600 11698 10238
40871 22537 110313 125531 110165 96729 83435 38881 58939 51639
Meshes A2 and F2 have been selectively refined near the spheres relative to meshes A1 and F1, respectively.
where IID,Y is the strain rate at yield, obtained by solving Eq. (6) at the yield stress: −1 2 1 IID,Y = 2 (−PB + PB−2 + 2BnP−1 B ) .
(9)
The scaled strain rate can be thought of as the relative state of yield; IID = 1 corresponds to the yield surfaces defined by the calculated stress for a fixed value of PB . As we increase the value of PB , the position of the yield surface changes; the true yield surfaces are found in the Bingham limit (PB → ∞). We note that IID increases in some regions while it decreases in others as we increase PB ; a representative profile is shown in Fig. 4. Moreover, this behavior is observed up to values of PB beyond which numerical results are unreliable. Extrapolating to the Bingham limit, we identify the regions where IID increases with PB as yielded, and the regions where IID tends to zero with increasing PB as unyielded. We can
Fig. 3. Contours of strain rate for two spheres falling collinearly calculated using (a) a coarse mesh (F2) and (b) a fine mesh (F1). L/R = 4.0, Bn = 340.7, PB = 5000. The contours range logarithmically from 10−14 to 106 , in 20 intervals.
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
55
Fig. 4. Profiles of the scaled strain rate (IID ) vs. radial position along the mid-plane (z = −L/2) for various values of PB . Two spheres falling collinearly: L/R = 4.0 and Bn = 340.7.
56
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
then define the extrapolated yield surfaces as the loci of points where IID is invariant with respect to PB . We can use this extrapolation to identify yielded regions even where the calculated scaled strain rates are less than unity (as in Fig. 4b); the sensitivity of the results to PB must be monitored, however, to ensure that the regions of increasing and decreasing IID are indeed tending to infinity and zero, respectively. For each value of L/R, simulations were run for different values of PB . We define a quantity φ≡
∂ log IID , ∂ log PB
(10)
which represents the change in relative state of yield with regularization parameter. The value of φ was calculated using a linear least squares fit at each nodal point and interpolated within the mesh. The extrapolated yield surfaces were defined by φ = 0; regions of positive and negative φ were identified as yielded and unyielded, respectively. Preliminary studies were performed using PB = 5 × 103 , 6 × 103 , 7 × 103 , 8 × 103 , 9 × 103 , and 104 ; the least squares fit using the linear approximation was very good (the coefficient of determination, r2 , was found to be greater than 0.95 over nearly the entire field). Moreover, the results using the two values PB = 5 × 103 and 8 × 103 were identical to the results using the full range of PB (5 × 103 , 6 × 103 , 7 × 103 , 8 × 103 , 9 × 103 , 104 ). We thus used PB = {5 × 103 , 8 × 103 } for the present calculations. Simulations using larger values of PB (up to 106 ) also showed good agreement, though numerical round-off error was readily apparent in some of these cases. The Stokes drag coefficient, defined CS =
Fdrag , 6πηp VR
(11)
where Fdrag is the drag experienced by each sphere, was also calculated for varying sphere separations. 4. Results: two spheres falling collinearly 4.1. Yield surfaces In calculating the locus of points where φ = 0, several different surfaces were identified, some of which may represent true extrapolated yield surfaces and some of which are mesh related artifacts. Fig. 5a shows the various surfaces for two spheres falling in a line with L/R = 7.0. As expected, we see an envelope surrounding each sphere (α), corresponding to the yield envelope observed in previous studies [2,4]. The yielded envelope is nearly identical to the envelope calculated for a single sphere in our previous work [1]. In addition, we see apparent yield surfaces near the corners of the mesh (β, γ) and also near the mid-plane between the spheres (δ). When we plot the values of the strain rate along the apparent yield surfaces, as in Fig. 5b, we find that the strain rates associated with surfaces β, γ, and δ, are several orders of magnitude lower than those on the yield envelope. Moreover, regions β and γ were dependent on the size of the computational domain (this behavior was also observed for single-sphere calculations), and the regions moved when the calculation was repeated on a larger mesh (Rmesh = Zmesh = 10). Region δ disappeared with selective mesh refinement. Examination of the strain rates near the mid-plane (Fig. 6) show that
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
57
Fig. 5. (a) Calculated yield surface for two spheres falling collinearly: Bn = 340.7 and L/R = 7.0, using mesh A1. Unyielded regions are shaded. β, γ, δ indicate anomalous yield surfaces. (b) Strain rate plotted against radial position along calculated yield surfaces; two spheres falling collinearly: Bn = 340.7 and L/R = 7.0.
58
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
Fig. 6. Profiles of the scaled strain rate (IID ) vs. radial position along the mid-plane (z = −L/2). Coarse mesh (A1) results in spurious yield surfaces. Two spheres falling collinearly: Bn = 340.7 and L/R = 7.0.
this apparent yielded region (δ) actually corresponds to the strain rate passing through zero; because this results in a cusp in IID , small mesh-dependent errors can cause the sign of φ to change, resulting in an apparent yield surface when the mesh is insufficiently refined. Hence, we clearly see that surfaces β, γ, and δ of Fig. 5a are mesh-related artifacts, and increased local mesh refinement can be used to eliminate these regions (the remaining meshes were refined to eliminate or minimize the mid-plane regions). We can also use the magnitude of the strain rate along the calculated yield surface to discriminate the true yield surfaces from these spurious mesh-related yield surfaces. This method is particularly useful for discriminating the corner regions, as it is computationally inefficient to refine the mesh regions far away from the spheres. Fig. 7 shows the true yield surfaces for two spheres falling in a line at smaller values of L/R (only the portion of the computational region near the yielded region is shown in this and subsequent figures). These surfaces were determined from the calculated yield surfaces by eliminating those surfaces where IID is less than some threshold value, as determined by a plot analogous to Fig. 5b. This threshold value is dependent on the local refinement of the mesh, and was found to be around 10−7 for all of the meshes except meshes G and H, where 10−5 was used. At L/R = 6.0 we see that the yield
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
59
Fig. 7. Calculated yield surfaces for two spheres falling collinearly: Bn = 340.7 and (a) L/R = 6.0; (b) L/R = 5.0; (c) L/R = 4.5; (d) L/R = 4.0; (e) L/R = 2.5; (f) L/R = 2.1. Spurious mesh-related yield surfaces have been removed.
envelopes are separate and essentially unchanged from those calculated for a single sphere. For these larger separations (L/R ≥ 6.0), there is essentially no interaction between the spheres. We can see that as we decrease the sphere separation, the spheres begin to interact and the yield envelopes eventually coalesce.
60
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
Fig. 8. Streamlines of the re-circulation pattern within the yielded envelope, two spheres falling collinearly: Bn = 340.7 and L/R = 4.0; 20 evenly spaced streamlines (ψ = −0.03–0.5). The axial velocity (vz ) profile along the mid-plane (z = −2.0) is shown for reference.
For L/R ≤ 5.0, the two spheres are no longer contained in separate yield envelopes, but rather in a single coalesced envelope. Although the yield surfaces for two isolated spheres would not intersect at L/R = 5, the stress field between the two spheres is not a linear superposition of two isolated spheres, so it is not surprising that the spheres are found to interact at this distance. We note distinct unyielded regions within the yielded envelope: along the mid-plane, centered at r ∼ 2; along the axis of symmetry, centered at the mid-plane between the spheres; and on the poles of the spheres. The region along the mid-plane corresponds to a plug-like flow of the material due to re-circulation within the yielded envelope. The re-circulation can be seen by examining the streamlines in Fig. 8 and is clearly required to preserve mass conservation within the yielded envelope. Note that the unyielded region is not in rotation, which would be inadmissible given the axisymmetric geometry, but rather moves uniformly opposite the motion of the spheres. This region must also exist as at least a point within the separate yield envelopes (L/R ≥ 6.0), although it may not be large enough to be visible in these calculations. The mid-plane unyielded region decreases in size with decreasing sphere separation and vanishes as the separation approaches 2R.
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
61
Fig. 9. Comparison of calculated yield surfaces for two spheres falling collinearly (solid lines) and a capped-cylinder (dashed lines): Bn = 340.7 and (a) L/R = 2.1; (b) L/R = 2.5; (c) L/R = 4.0; (d) L/R = 5.0.
The size of the coalesced envelope decreases with decreasing sphere separation while maintaining the same approximate shape and aspect ratio. By comparing the results with results for two spheres connected by a rigid bridge (i.e. a capped cylinder), as shown in Fig. 9, we see that the shape of the outer yield envelope is relatively insensitive to the shapes of the interior unyielded regions and appears to be dominated by the shape of the leading (and trailing) edge to the flow.
62
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
Fig. 10. Calculated yield surface for two spheres falling collinearly: Bn = 340.7 and L/R = 5.5. The annular region is expected to vanish in the Bingham limit.
The unyielded region along the axis of symmetry decreases in size with decreasing sphere separation, maintaining an approximately constant distance from the spheres, until the separation reaches about 4R; the region then becomes larger and connects the two spheres. The plug of unyielded material connecting the spheres then travels at the same velocity as the spheres. The polar caps probably exist on both the fore and aft surfaces of each sphere because of the stagnation points [2], but they are too small to be distinguishable in some of the simulations. The region on the pole closest to the mid-plane appears to increase in size with decreasing sphere separation until it abruptly coalesces when the spheres become connected at L/R = 4.0. The region on the opposite pole also appears to be increasing in size with decreasing sphere separation, but the regions are so small that the current calculations cannot resolve them accurately. The sphere separation of 5.5 (Fig. 10) represents a transition between separate yielded envelopes and a single, coalesced yielded envelope, in that the calculated yield surface appears as two separate envelopes connected by a thin annular channel. The axial velocities in this annular region and the region within it are in the same direction, however, and this configuration cannot conserve mass except through movement in the opposite direction in the surrounding unyielded region. Since the far-field velocities are fixed at zero, such a flow cannot occur in the Bingham limit, and the calculated velocities in the apparent plug do, in fact, decrease with increasing PB , as does the thickness of the annulus. A regularization parameter sufficiently large to see the complete vanishing of the annular channel was inaccessible for this transition case because of round-off error. The yield surfaces for smaller values of Bn were calculated for L/R = 4.0 and are shown in Fig. 11. The yield surfaces are relatively insensitive to Bn in this range; the height of the yielded envelope along the axis of symmetry (measured from the center of the nearest sphere) increases by only 14%, and the width at the mid-plane increases by only 11%, as Bn is decreased from 340.7 to 14.91. These results are similar to previous single-sphere calculations, which show slightly larger increases (20–30%) over the same range of Bn [1,2].
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
63
Fig. 11. Calculated yield surfaces for two spheres falling collinearly with L/R = 4.0 at various values of Bn.
4.2. Drag coefficients and variational integrals The calculated values of CS , H and K are insensitive to the value of PB ; they vary by less than 0.5% with a twofold change in the regularization parameter (5 × 103 –104 ). The accuracy of the results is confirmed by the small relative difference between H and K shown in Table 2. This is true even for the calculations at L/R = 5.5, where the anomalous annular channel connecting the yield surfaces appears (cf. Fig. 10), Table 2 Calculated variational integrals for two identical rigid spheres falling along the line of centers: Bn = 340.7 and PB = 104 L/R
Mesh
H
K
(H − K)/K (%)
7.0 7.0 6.0 5.5 5.0 4.5 4.0 4.0 2.5 2.1
A1 A2 B C D E F1 F2 G H
5020.7 5075.0 5032.3 5044.5 5014.5 4915.1 4729.4 4737.7 4037.5 3820.7
5016.9 4987.2 5000.1 5003.0 4975.7 4866.5 4635.9 4584.9 3864.1 3652.2
0.08 1.8 0.64 0.83 0.78 1.0 2.0 3.3 4.5 4.6
64
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67 1 0.9 Newtonian result
normalized drag coefficient (C S /CS,single sphere)
0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 2
3
4
5
6
7
8
9
10
sphere separation (L/R)
Fig. 12. Plot of normalized drag coefficient (CS /CS,single sphere ) for two spheres falling collinearly as a function of sphere separation, Bn = 340.7. Open symbols were calculated using less refined meshes.
and underscores the sensitivity of the location of the yield surfaces to the details of the strain field. The drag coefficients normalized by the drag coefficient for a single, isolated sphere are shown for various sphere separations in Fig. 12, together with the results for a Newtonian fluid [11,15]. For L/R ≥ 6, the drag experienced by each sphere is identical to that for an isolated sphere, as expected from the shapes of the yield surfaces (Figs. 7 and 8). The range of interaction in the yield stress material is much smaller than in the Newtonian case; the presence of the second sphere reduces the drag at small separations, but to a lesser extent than in a Newtonian fluid. An approximately 30% drag reduction is observed in the limiting case when the spheres are touching (L/R = 2.0). 5. Results: two approaching spheres Several sphere separations were examined for the case of two rigid spheres approaching at equal speed. The drag experienced by each sphere is tabulated in Table 3. There is no effect on the drag force for sphere separations greater than or equal to 4.5; this contrasts with the case of two spheres falling in the same direction, where interaction is observed up to L/R = 5. For smaller sphere separations, the spheres are increasingly affected by one another. Unlike the Newtonian fluid case, however, the spheres experience a slight drag reduction rather then an augmentation [12]; the reason for the net decrease in drag is that the interactions of the two spheres causes a decrease in the local viscosity because of the shear-thinning characteristic of the Bingham material.
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
65
Table 3 Calculated variational integrals and drag correction factor for two identical rigid spheres approaching one another at the same speed: Bn = 340.7 and PB = 5 × 103 L/R
Mesh
H
K
CS /CS,single sphere
4.5 4.0 2.5
E F1 G
5032.5 5016.1 4792.4
4999.7 4979.8 4736.6
1.00 0.99 0.94
Fig. 13. Calculated yield surfaces for two approaching spheres: Bn = 340.7 and (a) L/R = 4.5; (b) L/R = 4.0.
Fig. 14. Streamlines of the re-circulation pattern within the yielded envelope, two approaching spheres: Bn = 340.7 and L/R = 4.0; 20 evenly spaced streamlines (ψ = 0.01–0.5).
66
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
The extrapolated yield surfaces for two approaching spheres are shown in Fig. 13. At L/R = 4.5, the spheres are enclosed in separate yielded envelopes, as for a single, isolated sphere; this result is consistent with the lack of interaction observed in the drag calculations. At L/R = 4.0, the yielded envelopes are coalesced. The shape of the coalesced envelope is qualitatively different from that observed for two spheres falling in a line; these differences are highlighted by comparing the streamlines for each case (Fig. 14 for two approaching spheres and Fig. 8 for two spheres falling collinearly).
6. Conclusions As in our earlier study of a single sphere in a Bingham material [1], numerical error for calculations of two spheres in creeping flow along their line of centers precludes the use of regularization parameters sufficiently large to determine the location of the yield surfaces directly from the yield condition. The yield surfaces are therefore determined from the loci of values of the normalized second invariant of the deformation rate that are invariant with respect to the regularization parameter, together with refinement of the mesh and consideration of the magnitude of the deformation rate on the apparent yield surfaces. The minimum value of the regularization parameter needed to reflect ideal Bingham behavior is dependent on the details of the flow field. For sphere separations L/R ≥ 6, the yield surfaces are identical to those for single isolated spheres. The spheres interact at separations larger than those that would be predicted from single-sphere calculations; the yielded envelopes surrounding each sphere coalesce, and internal plug-like regions appear along the symmetry axis and along the mid-plane. A decrease in the drag coefficient of up to 30% accompanies the coalescence of the yielded envelopes. The size of the coalesced yield surface in insensitive to Bingham number. The range of interaction is shorter for two approaching spheres than for two spheres falling in a line, and the shapes of the coalesced yield surfaces are qualitatively different. Only a slight drag decrease compared to the single-sphere case is predicted, but this contrasts sharply with the drag increase that is seen in a Newtonian fluid.
Acknowledgements This research was supported in part by the Department of Energy, award number DE-AC03-76SF00098, and through computing resources provided by the National Partnership for Advanced Computational Infrastructure at the University of Michigan Center for Parallel Computing. Financial support through a National Science Foundation Graduate Research Fellowship to B.T.L. is also gratefully acknowledged.
References [1] B.T. Liu, S.J. Muller, M.M. Denn, Convergence of a regularization method for creeping flow of a Bingham material about a rigid sphere, J. Non-Newtonian Fluid Mech. 102 (2002) 179. [2] A.N. Beris, J.A. Tsamopoulos, R.C. Armstrong, R.A. Brown, Creeping motion of a sphere through a Bingham plastic, J. Fluid Mech. 158 (1985) 219. [3] N. Yoshioka, K. Adachi, H. Ishimura, On creeping flow of a viscoplastic fluid past a sphere, Kagaku Kogaku 35 (1971) 1144 (in Japanese).
B.T. Liu et al. / J. Non-Newtonian Fluid Mech. 113 (2003) 49–67
67
[4] J. Blackery, E. Mitsoulis, Creeping motion of a sphere in tubes filled with a Bingham plastic material, J. Non-Newtonian Fluid Mech. 70 (1997) 59. [5] E. Mitsoulis, Effect of rheological properties on the drag coefficient for creeping motion around a sphere falling in a tightly-fitting tube, J. Non-Newtonian Fluid Mech. 74 (1998) 263. [6] J. Oldroyd, A rational formulation of the equations of plastic flow for a Bingham solid, Proc. Camb. Philos. Soc. 43 (1947) 100. [7] J. Oldroyd, Two-dimensional plastic flow of a Bingham solid, Proc. Camb. Philos. Soc. 43 (1947) 383. [8] W. Prager, Introduction to Continuum Mechanics, Ginn, Boston, 1961. [9] R. Glowinski, J.L. Lions, R. Tremolieres, Numerical Analysis of Variational Inequalities, North-Holland, Amsterdam, 1981. [10] W. Prager, in: G. Birkhoff, G. Kuerti, G. Szegö (Eds.), Studies in Mathematics and Mechanics, Academic Press, New York, 1954, p. 208. [11] M. Stimson, G.B. Jeffery, The motion of two spheres in a viscous fluid, Proc. R. Soc. London A111 (1926) 110. [12] H. Brenner, The slow motion of a sphere through a viscous fluid towards a plane surface, Chem. Eng. Sci. 16 (1961) 242. [13] M. Bercovier, M. Engelman, A finite-element method for the numerical solution of viscous incompressible flows, J. Comp. Phys. 30 (1979) 181. [14] M. Bercovier, M. Engelman, A finite-element method for incompressible non-Newtonian flows, J. Comp. Phys. 36 (1980) 313. [15] J. Happel, R. Pfeffer, The motion of two spheres following each other in a viscous fluid, Am. Inst. Chem. Eng. J. 6 (1960) 129.