On the viscoplastic squeeze flow between two identical infinite circular cylinders
aa r X i v : . [ phy s i c s . f l u - dyn ] F e b On the viscoplastic squeeze flow between two identical infinitecircular cylinders
A.R. Koblitz, S. Lovett, and N. Nikiforakis Department of Physics, Cavendish Laboratory,J J Thomson Avenue, Cambridge, CB3 0HE, UK ∗ Schlumberger Gould Research Centre, High Cross,Madingley Road, Cambridge CB3 0EL, UK (Dated: February 8, 2018)
Abstract
Direct numerical simulations of closely interacting infinite circular cylinders in a Bingham fluidare presented, and results compared to asymptotic solutions based on lubrication theory in thegap. Unlike for a Newtonian fluid, the macroscopic flow outside of the gap between the cylindersis shown to have a large effect on the pressure profile within the gap and the resulting lubricationforce on the cylinders. The presented results indicate that the asymptotic lubrication solution canbe used to predict the lubrication pressure only if the surrounding viscoplastic matrix is yieldedby a macroscopic flow. This has implications for the use of sub-grid-scale lubrication models insimulations of non-colloidal particulate suspensions in viscoplastic fluids. ∗ A.R. Koblitz acknowledges financial support from the EPSRC Centre for Doctoral Training in Compu-tational Methods for Materials Science under grant EP/L015552/1; This work was supported by theSchlumberger Gould Research Centre.; [email protected] . INTRODUCTION Complex fluids are ubiquitous in natural and industrial processes, from food processing,to lava or debris flows, to oil and gas applications. The mechanical behaviour of thesefluids arises from the microstructure of the fluid, for example emulsion droplets and clays indrilling muds, or polymer chains in viscoelastic fluids. When non-colloidal particles muchlarger than the fluid microstructure are added, the system can be thought of as a particulatesuspension in a complex (continuum) fluid. Examples of these types of systems include freshconcrete and debris flows [1]. The hydrodynamic interaction between particles affects thesuspension bulk properties and dynamics and is of great interest. In the case of a Newtonianfluid, analytical solutions exist for slow flow past spheres and cylinders [2, 3] and the squeezeflow between them using asymptotic analysis. Viscoplastic fluids, of interest to this work,are characterised by a discontinuous nonlinear constitutive equation thereby introducingadditional complexities when analytical solutions are sought.So far, studies on interacting spheres and cylinders in viscoplastic flows have largelyfocussed on drag and pressure drop (in the case of flow past arrays) of collinear arrangements,aligned either parallel or perpendicular to the flow [4–10]. Numerical studies using theBingham constitutive law have been found to be in good agreement with experimental workusing Carbopol 940 gels, developing drag correlations and stability criteria (with respectto sedimentation) [4, 7, 9, 10]. Viscoplastic squeeze flow between coaxial cylindrical diskshas been studied analytically for both planar [11] and axisymmetric [12, 13] configurations.The configuration of collinearly approaching bodies in a viscoplastic flow has received onlycursory attention in numerical studies, eg Yu and Wachs [8], Tokpavi et al. [9], with noexamination of the interstitial squeeze flow.This study therefore examines the two-dimensional squeeze flow between two approachinginfinite circular cylinders in a Bingham viscoplastic fluid by direct numerical simulation. Theconfiguration studied is such that the gap between the two cylinders is small (1 % of thecylinder radius). We also make use of the asymptotic analysis by Balmforth [14] to computeleading order lubrication solutions for the squeeze flow between two approaching cylindersin a Bingham fluid. We compare the analytical and numerical solutions and demonstratethat in a quasi-unconfined system the squeeze flow is greatly affected by flow external tothe gap, but that the asymptotic solution may be recovered under certain flow conditionsin the wider domain. This is contrary to the Newtonian equivalent, and has implicationson using the viscoplastic lubrication force approximation as a sub-grid-scale model in coarsesimulation techniques.The paper is organised as follows. In section II we present the problem of interest andbriefly describe the solution strategy employed for the direct numerical simulations, and thelubrication theory calculations. In section III we present direct numerical simulations of thequasi-unconfined system. These are compared to simulations of the domain restricted tothe gap only, and to the asymptotic solutions from lubrication theory. These comparisonsdemonstrate the influence of the wider flow field on the lubrication pressure. In section IVwe discuss the results and the implications for sub-grid-scale modelling.
II. MATHEMATICAL FORMULATION AND SOLUTION
We consider the slow, steady flow of an incompressible viscoplastic fluid around two rigid,infinite circular cylinders. The fluid has velocity ˆ u (ˆ x ), pressure ˆ p (ˆ x ) and a symmetric total2 u = ( V / , u = ( − V / , D D H Ω Ω y x y x FIG. 1. Schematic showing the problem geometry for the entire flow system (Ω ) investigatedthrough numerical methods and the reduced system (Ω ) investigated with both analytical andnumerical methods. stress tensor ˆ τ − ˆ p δ , where variables with a hat are dimensional. In the absence of inertia,the conservation of mass is ∂ ˆ u∂ ˆ x + ∂ ˆ v∂ ˆ y = 0 , (1)and the conservation of momentum is ∂ ˆ τ xx ∂ ˆ x + ∂ ˆ τ xy ∂ ˆ y − ∂ ˆ p∂ ˆ x = 0 , (2) ∂ ˆ τ yx ∂ ˆ x + ∂ ˆ τ yy ∂ ˆ y − ∂ ˆ p∂ ˆ y = 0 . (3)As a constitutive law we use the Bingham model ( ˆ τ ij = (cid:16) η + ˆ τ Y ˆ˙ γ (cid:17) ˆ˙ γ ij if ˆ τ > ˆ τ Y ,ˆ˙ γ ij = 0 if ˆ τ ≤ ˆ τ Y , (4)where ˆ τ Y and ˆ η are the yield stress and the plastic viscosity of the fluid, respectively, ˆ˙ γ ij isthe rate of strain tensor associated with the velocity field, andˆ˙ γ ij = 12 (cid:18) ∂ ˆ u i ∂ ˆ x j + ∂ ˆ u j ∂ ˆ x i (cid:19) , ˆ˙ γ = r
12 ˆ˙ γ ij ˆ˙ γ ij , ˆ τ = r
12 ˆ τ ij ˆ τ ij . (5)The problem geometries are depicted in figure 1, where the inset highlights the portionof the system considered in the analytical investigation. Aligning the system mid-plane ina Cartesian coordinate system the two cylinders are placed with their centres located at( − H/ − D/ ,
0) and ( H/ D/ , H is the minimum separation distance and D the cylinder diameter. The computational domain for the whole system has dimensions10 D × D , which is sufficiently large for the cylinders to be essentially unconfined: waningstresses away from the moving cylinders lead to the formation of a yield envelope in theimmediate vicinity of the cylinders, outside of which the fluid forms a rigid plug attachedto the domain walls. Because the fluid in the far field is unyielded for the range of yield3tresses explored in this study, we set no-slip boundaries at y ± . D and pressure inletsand outlets at x ± D . The cylinders have a constant relative approach velocity of V . A. Large-scale non-dimensionalisation
Choosing a velocity scale of V , length scale of D , shear rate scale of V /D , and stressscale of ˆ ηV /D , we obtain the dimensionless equations ∂u∂x + ∂v∂y = 0 , (6) ∂τ xx ∂x + ∂τ xy ∂y − ∂p∂x = 0 , (7) ∂τ yx ∂x + ∂τ yy ∂y − ∂p∂y = 0 , (8) ( τ ij = (cid:16) Bn ˙ γ (cid:17) ˙ γ ij if τ > Bn ,˙ γ ij = 0 if τ ≤ Bn , (9)where Bn := ˆ τ Y D ˆ ηV (10)is a Bingham number for the macroscopic flow external to the gap. B. Computational method
We numerically compute the solution of equations 6–9 for two approaching cylinders witha small gap size (
H/R = 0 .
01 where R ≡ D/ et al. [17] where its efficacy for particulate flow simulationswas demonstrated. Briefly, the overset grid method represents a complex domain usingmultiple body-fitted curvilinear grids that are allowed to overlap whilst being logically rect-angular. The overlapping aspect brings flexibility and efficiency to grid generation, which isbeneficial for moving body problems. Here, since the cylinders are static, the chief benefit ofthe overset grid method is that the grids can be locally refined near the gap whilst keepingthe grids logically rectangular. The resultant linear systems are solved using the MUMPSlibrary [18], a massively parallel direct linear solver. We use meshes with a minimum of 15points across the narrowest part of the gap and cluster grid points near the cylinder surfacesand wider gap region by stretching the constituent grids.Applying a standard finite difference method to equations (6)–(9) is not straightforward,due to the non-differentiable plastic dissipation term. A straight-forward way of dealing withthis numerical difficulty is to regularize equation (9) by removing the singularity at ˙ γ = 0.This approach has been used in studies of viscoplastic flows past bluff bodies, see [19–21].4owever this can yield inaccurate results, especially for lubrication-type flows or if flowstability or finite-time stoppage are of critical interest [22–24]. Instead, we use an iterativemethod based on the variational form of the Bingham problem, established by Duvaut andLions [25], which forms the basis for the widely used augmented Lagrangian (AL) firstproposed by Glowinski [26]. This formulation is commonly known as ALG2 and is usedextensively in the literature, see Yu and Wachs [8], Muravleva [11], Chaparian and Frigaard[27] and references therein, so we do not give details here. For its solution we use the Uzawatype algorithm of Olshanskii [28] and Muravleva and Olshanskii [29]. C. Lubrication flow in the gap
The problem shown in the inset of figure 1, i.e. the narrow gap between two symmetricsurfaces approaching with relative speed V , has an asymptotic solution due to Balmforth[14], if the gap H is small compared to the cylinder radius R . In this section we give anoverview of this solution; in section III we will compare this to fully numerical solutionsboth in the restricted domain (inset of figure 1) and the full domain. Note that this sectionconsiders a non-dimensionalisation of the governing equations appropriate to the gap scale;the non-dimensionalisation given previously in section II A is appropriate for the macroscopicflow. We take x to be the coordinate across the gap and y the coordinate along the gap,consistent with the setup shown in figure 1.We write ˆ u ≡ (ˆ u, ˆ v ) and without loss of generalityˆ τ ≡ (cid:18) ˆ σ ˆ ψ ˆ ψ − ˆ σ (cid:19) . (11)Following the approach in [14], variables are scaled as x = ˆ x/ H , y = ˆ y/ L , u = ˆ u/ U , v = ˆ v/ ( U /ǫ ) , p = ˆ p/ P , (12)where ǫ ≡ H / L is a small parameter. This implies the scaled continuity equation is ∂u∂x + ∂v∂y = 0 . (13)The stress scale is chosen as τ = ˆ τ / ( ǫ P ), which implies ∂p∂x = ǫ ∂σ∂x + ǫ ∂ψ∂y , ∂p∂y = ∂ψ∂x + ǫ ∂σ∂y , (14)so that the main force balance (to O ( ǫ )) is between the axial pressure gradient and transverseshear stress gradient. Strain rates are scaled by ( U /ǫ ) / H , giving˙ γ = s (cid:18) ǫ ∂u∂y + ∂v∂x (cid:19) + ǫ (cid:18) ∂u∂x (cid:19) . (15)5he above scaling implies in the yielded regions τ ij = (cid:18) η U ǫ PH + ˆ τ Y ǫ P ˙ γ (cid:19) ˙ γ ij . (16)The velocity scale is set by the motion of the cylinders as U := V and therefore the pressurescale is chosen as P := ˆ ηVǫ H . (17)We additionally fix the characteristic length and gap scales as L = R and H = H , respec-tively. This gives the scaled constitutive equation as ( τ ij = (cid:16) B ∗ ˙ γ (cid:17) ˙ γ ij if τ > B ∗ ,˙ γ ij = 0 if τ ≤ B ∗ , (18)where B ∗ := ˆ τ Y ǫ P (19)is a Bingham number for the squeeze flow in the gap. Note that B ∗ / Bn = ǫ /
2; the squeezeflow ‘sees’ a much lower Bingham number than the macroscopic flow around the cylinders.
1. Leading-order solution
The components of the shear rate tensor are ψ ≡ ˙ γ xy = 12 (cid:18) ǫ ∂u∂y + ∂v∂x (cid:19) , (20) σ ≡ ˙ γ xx = ǫ ∂u∂x . (21)Therefore, discarding terms of O ( ǫ ), the shear rate magnitude is˙ γ = 12 (cid:12)(cid:12)(cid:12)(cid:12) ∂v∂x (cid:12)(cid:12)(cid:12)(cid:12) , (22)and in the fully yielded part of the flow ψ ≫ σ . Equation 18 is used to write ψ = ∂v∂x + B ∗ sgn (cid:18) ∂v∂x (cid:19) , (23)and the main force balance reduces to ∂p∂x = 0 ⇒ p = p ( y ) , ∂p∂y = ∂ψ∂x ⇒ ψ = x ∂p∂y , (24)the constant vanishing by symmetry, meaning that the pressure gradient is, to leading order,constant across the gap and balanced along the gap by the transverse shear stress. Exploitingthe symmetry of the configuration, in the quadrant x > y > v > v∂x <
0, and so from the main force balance and constitutive law we find the velocity profileacross the gap ∂v∂x = x ∂p∂y + B ∗ , (25)which may be integrated to give v = ( − ∂p∂y (cid:0) h − x (cid:1) (cid:0) h − X + x (cid:1) , X < x ≤ h ( y ) − ∂p∂y (cid:0) h − X (cid:1) , ≤ x ≤ X, (26)where X ≡ B ∗ / | ∂p∂y | is the plug boundary location, and we have used a no-slip boundary con-dition at the cylinder surface, located at x = h ( y ). The continuity equation and boundaryconditions imply a flow rate constraint ∂∂y Z h − h v d x = 1 (27)which, when evaluated using the velocity solution, gives a cubic equation for the pressuregradient ∂p∂y ( y ): − ∂p∂y ( h + X )( h − X ) = y. (28)It can be shown that the plug in the region | x | < X undergoes O ( ǫ ) plastic flow, which isnot present in the above asymptotic solution. This may be recovered by keeping terms O ( ǫ )and is sometimes referred to as a pseudo-plug; it does not change the equation for thepressure gradient to leading order [14].For two converging cylinders the non-dimensional separation distance is h ( y ) = 1 + 2 ǫ (cid:16) − p − y (cid:17) , ≤ | y | < . (29)We numerically evaluate equation 28 to compute ∂p∂y ( y ) and thence p ( y ), with an additionalambient pressure constraint outside the disks enforced as p (1) = 0. The leading-orderlubrication force is then numerically computed as 2 R p d y . D. Flow field diagnostics
In order to classify the structure of the numerically-calculated flowfields we make use ofan invariant measure of the velocity gradient tensor that gives an indication of the relativestrength of the shear rate tensor and vorticity field [30] Q = − ∂ ˆ u i ∂ ˆ x j ∂ ˆ u j ∂ ˆ x i = − (cid:18) ˆ˙ γ ij ˆ˙ γ ij −
12 ˆ ω (cid:19) , (30)7here ˆ ω is the vorticity. We use the normalized form of (30)Λ = ˆ˙ γ ij ˆ˙ γ ij − /
2( ˆ ω )ˆ˙ γ ij ˆ˙ γ ij + 1 /
2( ˆ ω ) , (31)such that values of Λ = − , , W , is calculated by integrating the rate of viscous dissipa-tion, ˆΦ = ˆ τ ij ˆ˙ γ ij , over a suitable control volumeˆ˙ W (Ω) = Z Ω − V C ˆ τ ij ˆ˙ γ ij d V, (32)where V C is the volume occupied by the cylinders. This is scaled by the force on thecylinders and the closing velocity, W = F V , while the viscous dissipation is scaled using acharacteristic energy density scale E = ηV / H . III. RESULTS
We investigate the squeeze flow between two infinite circular cylinders in three differentcases based on the set-up shown in figure 1. Non-dimensionalisation is as described insection II A. The external Bingham number Bn is varied between 0 and 2000 in all cases,with the minimum separation distance kept constant at 0 .
01 non-dimensional units (i.e. 1%of the cylinder radius), resulting in a gap Bingham number B ∗ range of 0 to 0.1.Two cases use the full computational domain, labelled as Ω in figure 1, with differing farfield boundary conditions. The quiescent case (meaning here that the flow is zero outsidea finite yield envelope) has velocity outlets at the vertical domain boundaries, imposingan ambient pressure of p = 0. No-slip and no-penetration conditions are imposed on thehorizontal domain boundaries, allowing the cylinders to be surrounded by a bounded yieldedregion enclosed by a yield envelope.The shear flow case considers the same geometry as the quiescent case but with theintroduction of a macroscopic flow to raise the stress above the material yield stress in thefar field, thus removing the yield envelope. This is done by imposing wall velocities ± U w onthe horizontal domain boundaries, resulting in a macroscopic shear rate of ˙ γ = 5.Finally, we consider the reduced domain labelled as Ω in figure 1, only including the gapbetween the cylinders. No-slip and no-penetration conditions are applied on the cylindersurfaces, and symmetry conditions at x = ± IG. 2. Velocity magnitude plots for increasing Bn (left to right) with quiescent conditions (toprow) and macroscopic shear rate ˙ γ = 5 (bottom row).FIG. 3. Pressure contours for the two cylinder system in macroscopic shear flow (left) and quiescent(right) with Bn = 1000. A. Flow field kinematics
Figure 2 shows a binary yielded/unyielded mask in grey overlaid on colour maps ofthe velocity magnitude for the quiescent and sheared systems, with the Bingham numberincreasing from left to right. The unyielded regions are identified as areas where the secondinvariant of the shear stress falls below the yield stress, plus some small constant which wetake as 0 . Bn ≥
50 classicalfeatures of moving bodies in yield stress fluids can be seen: unyielded caps on the stagnation9
IG. 4. Left: Ratio of centreline pressure to surface pressure as a function of the local gap widthfor Bn = 0 (circles), 50 (diamonds), 500 (crosses), 1000 (triangles), and 2000 (squares). Right:Pressure distributions over the entire cylinder surface for the limiting Bn = 0 and Bn = 2000quiescent flow cases as a function of angle away from the gap centre. points, unyielded plugs in the equatorial planes of the cylinders, and a yield envelope fullysurrounding the two cylinder system [19, 27, 33–36]. As the Bingham number increases theunyielded stagnation caps and the equatorial plugs grow while the yield envelope shrinks.The bottom row of figure 2 corresponds to the shear flow case, where the backgroundshear flow has noticeably changed the yield surface features: stagnation points have shifted,leading to two caps on the rear of the cylinders, placed symmetrically about the longitudinalaxis, and one on the front of each cylinder. The equatorial plugs are no longer present, buttwo unattached plugs have formed in the gap openings, placed asymmetrically about thelongitudinal axis. For Bn >
50 central plugs can be seen fore and aft of the two cylindersystem.Figure 3 show contour plots of the pressure field for a quiescent (right panel) and shearflow (left panel) case at Bn = 1000, with the contours drawn at the same levels in bothpanels. The quiescent case shows a pressure drop from the gap to the rear stagnation cap,with roughly equally spaced iso-contours along shear layers attached to the cylinder surfacesand along the yield envelope boundary. In contrast, the shear flow case shows a rapidpressure decay along the gap with a more uniform pressure field outside of the gap. B. Small gap approximation
In section II C a lubrication approximation was constructed which has, to leading order,pressure constant across the gap. This approximation relies on the gap being small, viz. ǫ ≪
1, which is satisfied at the gap centre. However, since the approaching surfaces are elliptic(see equation (29)), this condition will be violated towards the gap exit.We investigate the validity of this constant pressure solution by examining the ratio ofthe pressure at the gap centre to that at the surface of the cylinder as a function of the localgap width. In the left panel of figure 4 we plot this gap-to-surface pressure ratio against the10ormalized gap width as a function of y as markers for Bn = 0–2000 ( B ∗ = 0–0 . . h (ˆ y ) / L > .
1. In the Newtonian case, the pressure ratiothen rapidly decreases away from the gap centre, becoming negligible at the gap exit. Allviscoplastic cases show similar behaviour to one another (as indicated by the marker overlapin the left panel of figure 4). The gap-to-surface pressure ratio remains close to unity closeto the gap centre before slowly decreasing. Unlike for the Newtonian case, no rapid decreaseis found when ˆ h (ˆ y ) / L > . . θ = 0 and the gap exits at θ = ( − π/ , π/ . C. Pressure profiles in the gap
Figure 5 presents pressure profiles through the centre of the gap, i.e. along the axis ofsymmetry. The direct numerical simulation (DNS) in the reduced domain gives a pressureprofile in excellent agreement with the DNS of the full system in the macroscopic shearflow. Both these cases are in good agreement with the asymptotic solution from lubricationtheory: peak pressures in the centre of the gap match well for the full Bn range explored.At higher Bn , the DNS pressures of the shear flow and reduced domain cases remain inagreement but decay more slowly than the asymptotic solution as the gap widens up; thisis where the lubrication approximation no longer holds.The pressure profiles for the DNS of the full system in the quiescent case are markedlydifferent to the asymptotic solution for Bn >
50: higher peak pressures and slower pressuredecay are evident, as is an exit pressure significantly higher than the ambient pressure (whichis 0).The relative change in peak pressure and pressure decay for increasing Bn discussedabove is evident in the surface pressure distribution shown in figure 4. Moreover, it isevident that the pressure contribution outside of the nominal gap region is negligible. Notethat for Bn = 2000 the pressure profiles look somewhat similar in magnitude between thequiescent and sheared cases. In fact their integrals differ by about a factor of two, implyinga factor of two difference in the repulsive force; this is discussed next.Figure 6 presents stacked area plots of the total drag force exerted on a cylinder, decom-posed into pressure and viscous contributions. From the left panel it is clear that the force11 IG. 5. Pressure profile along the gap centreline for quiescent (diamonds), sheared (circles), andreduced domain (square) systems with viscoplastic lubrication solution overlayed (solid line) with Bn = 50 , , , a ) through ( d ), respectively. on the cylinders in quiescent fluid is dramatically underestimated when the lubrication flowin the gap is considered in isolation. Both viscous and pressure contributions increase with Bn , but the viscous contribution remains small compared to that of the pressure. Discount-ing the viscous friction, the pressure alone—which remains localised to the gap—causes amore than two-fold increase in the drag force over the predictions from lubrication theory.However when a macroscopic shear flow is added (the right panel in figure 6), the totaldrag force is close to the asymptotic solution. This mirrors the trend found in the pressureprofiles in figure 5.Figure 7 shows local shear rates, ˙ γ local , for both the sheared and quiescent cases for thefull range of Bn . We define the local shear rate as an average over a H × H area inthe centre of the gap. The dramatic increase in pressure and drag force is not reflectedin the local shear rate: the local shear rates in the quiescent and sheared cases remain inclose agreement throughout the Bn range explored. From this we can conclude that themacroscopic flow does not affect the velocity field in the gap. The above computations werealso performed with a macroscopic shear rate one order of magnitude higher, showing noappreciable differences in the pressure and force results discussed above.12 IG. 6. Stacked area plots of the total drag force on a cylinder as a function of Bn for ( a ) quiescentand ( b ) shear flow background conditions. Overlaid are the predictions from viscoplastic lubricationtheory (squares).FIG. 7. Local shear rates in the gap centre for quiescent (squares) and sheared (diamonds) systemsfor Bn = 0 , , , , IG. 8. (a) Colour map of log ( ˙ γ ) and (b) normalized second invariants of the velocity gradienttensor with unyielded areas masked in grey with Bn = 1000 (quiescent case).FIG. 9. Left: rate of mechanical dissipation per unit Λ. Right: stacked area plot showing the rate ofwork done by viscous dissipation in the shear and plastic regions, scaled by W = F V , with markersindicating the power required to move the cylinders with the approach velocity (diamonds), thetotal viscous dissipation in the system (dashed line), the viscous dissipation in the gap region ofthe full system (squares) and the total viscous dissipation in the reduced system (triangles).
D. Viscous dissipation
The left panel of 8 shows a colour map of log ( ˙ γ ) in a quiescent case, while the rightpanel shows a contour plot of the normalized second invariant of the velocity gradient tensor,indicating where the fluid is irrotational (red), rotational (blue) and being sheared (green),14ith the yield surface overlaid. As expected, the plugs either side of the cylinders undergorigid body rotation. Strain rates are highest in the thin shear layers along the cylindersurface, along the yield envelope wall, and surrounding the jet of fluid squeezed out of thegap. Strain dominated regions are found in the core of the fluid jet squeezed out of the gap,and in the regions between the rotating plugs and the yield envelope. The strain rate inthese plastic flow regions is orders of magnitude lower than in the adjacent shear layers.We turn now to the energy dissipation in the fluid. The left panel of figure 9 showsthe rate of mechanical dissipation as a function of the topology parameter Λ. As would beexpected, no dissipation occurs in the regions undergoing rigid body rotation (Λ = − Bn > − / ≤ Λ ≤ / > /
3, respectively [31]. While the dissipation inplastic flow is small compared to that in shear flow, it still forms a significant source ofviscous dissipation outside of the gap area due to the size of the plastic flow regions (seefigure 8). Finally, the rate of work done in the gap region is very similar for both the full(Ω ) and reduced (Ω ) domains. This shows that the excess drag force on the cylinders inthe full domain compared to lubrication theory (figure 6( a )) arises from energy dissipationexternal to the gap. IV. CONCLUSIONS
In this paper we have presented results on the squeeze flow between two infinite circularcylinders in a Bingham fluid, which we use as a simple model for the flow of non-colloidalparticles in a viscoplastic fluid. Understanding this flow is essential to building models ofthe large-scale flow of such suspensions. Although the calculations presented here have beentwo-dimensional, we expect similar phenomena will occur in three dimensions (where theparticles would be spheres).In section III we presented results from three numerical experiments: two modelling theapproaching cylinders within a quiescent and a sheared fluid, and one modelling just thegap between the approaching cylinders, removing any external influence. We showed thatunlike for a Newtonian fluid, the macroscopic flow external to the gap has a large effect onthe lubrication forces felt by two cylinders in near-contact. In a quiescent Bingham fluid,the lubrication forces were approximately double those predicted by viscoplastic lubricationtheory, but were still caused primarily by the localised high lubrication pressure in thegap, as for a Newtonian fluid. The high lubrication pressure compared to theory is dueto the enclosing yield envelope which forms around the two particle system and causes arecirculating flow, introducing significant viscous dissipation into the system. Most of theextra viscous dissipation occurs in shear layers along the cylinder surface and yield envelopewalls, although at high Bingham numbers the contribution from plastic flow regions nearthe yield envelope becomes appreciable. 15ntroducing a macroscopic shear flow or modelling just the gap area between the cylindersgave nearly identical results, and agreed closely with the predictions from lubrication theory.We conclude that the background shear flow acts to eliminate the yield envelope in themacroscopic flow around the particles. This in turn removes the recirculating flow andcomplex flow structures, where large sources of viscous dissipation in the quiescent caseappear, and hence lowers the lubrication pressure and resulting lubrication force. Theresulting pressure profiles in the gap are well-described by lubrication theory local to thegap. The results indicate that the macroscopic shear rate does not appreciably affect thevelocity field in the narrow gap region. This conclusion is insensitive to the exact macroscopicshear rate used, provided the yield envelope is removed.The above implies that lubrication force models using an effective viscosity based on thelocal shear rate (such as the approach used for shear-thinning fluids in V´azquez-Quesada et al. [38]) may not be accurate for viscoplastic fluids. Instead, we suggest the use ofsub-grid-scale lubrication force models based on viscoplastic lubrication theory, with theunderstanding that they may become invalid in regions without a macroscopic stress abovethe yield stress, i.e. where particles become confined by their own yield envelopes. Thiswill allow for a large range of validity, for example in simulations of the type consideredin [38, 39] among others, where a dense suspension is subject to shear, and a sub-grid-scalelubrication force model is needed due to the close particle-particle approaches. Howeverin other cases, for example dilute particulate suspensions sedimenting in a quiescent fluid,we have shown in this paper that a sub-grid-scale lubrication force model based solely onlubrication theory in the gap may not be appropriate. Until a more sophisticated sub-grid-scale model is developed, the only current option is DNS computations with sufficiently highresolution in the inter-particle gaps. [1] G. Ovarlez, F. Mahaut, S. Deboeuf, N. Lenoir, S. Hormozi, and X. Chateau, “Flows ofsuspensions of particles in yield stress fluids,” J.Rheol. , 1449–1486 (2015).[2] M. Stimson and G. B. Jeffery, “The motion of two spheres in a viscous fluid,” P.Roy.Soc.A-Math.Phy. , 110–116 (1926).[3] A. Umemura, “Matched-asymptotic analysis of low-Reynolds-number flow past two equal cir-cular cylinders,” J.Fluid Mech. , 345–363 (1982).[4] B. T. Liu, S. J. Muller, and M. M. Denn, “Interactions of two rigid spheres translatingcollinearly in creeping flow in a Bingham material,” J.Non-Newtonian Fluid Mech. , 49–67 (2003).[5] M. R. Horsley, R. R. Horsley, K. C. Wilson, and R. L. Jones, “Non-Newtonian effects on fallvelocities of pairs of vertically aligned spheres,” J.Non-Newtonian Fluid Mech. , 147–152(2004).[6] P. Jie and Z. Ke-Qin, “Drag force of interacting coaxial spheres in viscoplastic fluids,” J.Non-Newtonian Fluid Mech. , 83–91 (2006).[7] O. Merkak, L. Jossic, and A. Magnin, “Spheres and interactions between spheres moving atvery low velocities in a yield stress fluid,” J.Non-Newtonian Fluid Mech. , 99–108 (2006).[8] Z. Yu and A. Wachs, “A fictitious domain method for dynamic simulation of particle sedi-mentation in Bingham fluids,” J.Non-Newtonian Fluid Mech. , 78–91 (2007).[9] D. L. Tokpavi, P. Jay, and A. Magnin, “Interaction between two circular cylinders in slow ow of Bingham viscoplastic fluid,” J.Non-Newtonian Fluid Mech. , 175–187 (2009).[10] L. Jossic and A. Magnin, “Drag of an isolated cylinder and interactions between two cylindersin yield stress fluids,” J.Non-Newtonian Fluid Mech. , 9–16 (2009).[11] L. Muravleva, “Squeeze plane flow of viscoplastic Bingham material,” J.Non-Newtonian FluidMech. , 148–161 (2015).[12] D. N. Smyrnaios and J. A. Tsamopoulos, “Squeeze flow of Bingham plastics,” J.Non-Newtonian Fluid Mech. , 165–190 (2001).[13] L. Muravleva, “Axisymmetric squeeze flow of a viscoplastic Bingham medium,”J.Non-Newtonian Fluid Mech. (2017), 10.1016/j.jnnfm.2017.09.006.[14] N. J. Balmforth, “Viscoplastic asymptotics and other techniques,” in Viscoplastic Fluids:From Theory to Application (Springer, 2017).[15] G. S. Chesshire and W. D. Henshaw, “Composite overlapping meshes for the solution of partialdifferential equations,” J.Comp.Phys. , 1–64 (1990).[16] W. D. Henshaw, Ogen: an overlapping grid generator for Overture , Research Report UCRL-MA-132237 (Lawrence Livermore National Laboratory, 1998).[17] A. R. Koblitz, S. Lovett, N. Nikiforakis, and W. D. Henshaw, “Direct numerical simulationof particulate flows with an overset grid method,” J.Comp.Phys. , 414–431 (2017).[18] P. R. Amestoy, I. S. Duff, J. L´Excellent, and J. Koster, “A fully asynchronous multifrontalsolver using distributed dynamic scheduling,” SIAM J.Matrix Anal.Appl. , 15–41 (2001).[19] D. L. Tokpavi, A. Magnin, and P. Jay, “Very slow flow of Bingham viscoplastic fluid arounda circular cylinder,” J.Non-Newtonian Fluid Mech. , 65–76 (2008).[20] T. Zisis and E. Mitsoulis, “Viscoplastic flow around a cylinder kept between parallel plates,”J.Non-Newtonian Fluid Mech. , 1–20 (2002).[21] S. Mossaz, P. Jay, and A. Magnin, “Criteria for the appearance of recirculating and non-stationary regimes behind in a viscoplastic fluid,” J.Non-Newtonian Fluid Mech. , 1525–1535 (2010).[22] I. A. Frigaard and C. Nouar, “On the usage of viscosity regularisation methods for visco-plasticfluid flow computation,” J.Non-Newtonian Fluid Mech. , 1–26 (2005).[23] A. Putz, I. A. Frigaard, and D. M. Martinez, “On the lubrication paradox and the useof regularisation methods for lubrication flows,” J.Non-Newtonian Fluid Mech. , 62–77(2009).[24] A. Wachs and I. A. Frigaard, “Particle settling in yield stress fluids: limiting time, distanceand applications,” J.Non-Newtonian Fluid Mech. , 189–204 (2016).[25] G. Duvaut and J. L. Lions, Les in´equations en m´ecanique et en physique (Dunod, 1972).[26] R. Glowinski,
Numerical methods for non-linear variational problems (Springer Verlag, 1984).[27] E. Chaparian and I. A. Frigaard, “Yield limit analysis of particles motion in a yield-stressfluid,” J.Fluid Mech. , 311–351 (2017).[28] M A. Olshanskii, “Analysis of semi-staggered finite-difference method with application toBingham flows,” Comput.Methods Appl.Mech.Engrg. , 975–985 (2009).[29] E. A. Muravleva and M. A. Olshanskii, “Two finite-difference schemes for calculation of Bing-ham fluid flows in a cavity,” Russ.J.Numer.Anal.Math.Modelling , 615–634 (2008).[30] P. A. Davidson, Turbulence: An Introduction for Scientists and Engineers (Oxford UniversityPress, 2004).[31] S. De, J.A. M. Kuipers, E.A.J. F. Peters, and J. T. Padding, “Viscoelastic flow simulationsin model porous media,” Phys.Rev.Fluids (2017).[32] I. A. Frigaard and D. P. Ryan, “Flow of a visco-plastic fluid in a channel of slowly varying idth,” J.Non-Newonian Fluid Mech. , 67–83 (2004).[33] A. Putz and I. A. Frigaard, “Creeping flow around particles in a Bingham fluid,” J.Non-Newtonian Fluid Mech. , 263–280 (2010).[34] A. N. Beris, J. H. Tsamopoulos, R. C. Armstrong, and R. A. Brown, “Creeping motion of asphere through a Bingham plastic,” J.Fluid Mech. , 219–244 (1985).[35] R. W. Ansley and T. N. Smith, “Motion of spherical particles in a Bingham plastic,” AIChE , 1193–1196 (1967).[36] K. Adachi and N. Yoshioka, “On creeping flow of a visco-plastic fluid past a circular cylinder,”Chem.Eng.Sci. , 215–226 (1973).[37] I. C. Walton and S. H. Bittleston, “The axial flow of a Bingham plastic in a narrow eccentricannulus,” J.Fluid Mech. , 39–60 (1991).[38] A. V´azquez-Quesada, R. I. Tanner, and M. Ellero, “Shear thinning of noncolloidal suspen-sions,” Phys.Rev.Lett. (2016).[39] X. Bian and M. Ellero, “A splitting integration scheme for the SPH simulation of concentratedparticle suspensions,” Comput.Phys.Comm. , 53–62 (2014)., 53–62 (2014).