On the Physics of Viscoplastic Fluid Flow in Non-circular Tubes

International Journal of Non–Linear Mechanics 88 (2017) 1–10 Contents lists available at ScienceDirect International Journal of Non–Linear Mechanics journal homepage: www.elsevier.com/locate/nlm On the physics of viscoplastic fluid flow in non-circular tubes a Mario F. Letelier , Dennis A. Siginer a,b,⁎ a crossmark , Cristian Barrera Hinojosa a Centro de Investigación en Creatividad y Educación Superior & Departmento de Ingeniería Mechánica, Universidad de Santiago de Chile, Santiago, Chile b Department of Mathematics and Statistical Sciences & Department of Mechanical, Energy and Industrial Engineering, Bostwana International University of Science and Technology, Palapye, Bostwana A R T I C L E I N F O A BS T RAC T Keywords: Viscoplastic Bingham Hershey-Bulkley Non-circular tube Plug and stagnant zones Flow of Bingham plastics through straight, long tubes is studied by means of a versatile analytical method that allows extending the study to a large range of tube geometries. The equation of motion is solved for general noncircular cross-sections obtained via a continuous and one-to-one mapping called the shape factor method. In particular the velocity field and associated plug and stagnant zones in tubes with equilateral triangular and square cross-section are explored. Shear stress normal to equal velocity lines, energy dissipation distribution and rate of flow are determined. Shear-thinning and shear-thickening effects on the flow, which cannot be accounted for with the Bingham model, are investigated using the Hershey-Bulkley constitutive formulation an extension of the Bingham model. The existence and the extent of undeformed regions in the flow field in a tube with equilateral triangular cross-section are predicted in the presence of shear-thinning and shear-thickening as a specific example. The mathematical flexibility of the analytical method allows the formulation of general results related to viscoplastic fluid flow with implications related to the design and optimization of physical systems for viscoplastic material transport and processing. 1. Introduction Knowledge of the flow of viscoplastic materials is relevant in many contexts such as flow of paints, pastes, suspensions in complex geometries with industrial applications, coating and mining, foodstuffs processing, cosmetic, and pharmaceutical and construction industries, ceramics extrusion, blood and other biological fluid flows, semi-solid materials and in some natural flows such as mud, lava displacements and debris flow. In all these applications as well as natural phenomena the rate of flow, the velocity distribution and the energy dissipation are important flow variables to determine. One model of viscoplastic fluid widely used is Bingham model for its capacity of predicting useful results in most areas of interest. The Bingham model becomes nonlinear for flow configurations different from parallel axisymmetric or unbounded parallel surfaces and, moreover, requires careful interpretation and analysis of mathematical results, which are meaningful only when all physically relevant conditions between stress and rate of deformation are met. Other types of non- Newtonian bounded flows, such as viscoelastic fluid flows in tubes, may be fully described physically over the whole flow region by means of mathematical or numerical results derived from the constitutive and linear momentum balance equations. But this is not the case for viscoplastic flows described by the Bingham model since it predicts physically mean- ⁎ ingless results in some zones that must be identified and characterized as plug zones and stagnant zones where there is no deformation. This is not explicitly predicted by the Bingham model, and must be deduced from conditions associated with the yield stress, tube contour and the related physical considerations. Understanding the dynamics of the formation of dead regions for instance is important to the design of extrusion geometries. It is quite difficult to model viscoplastic fluid flow and design operating systems in most real-life contexts. In particular the determination of the location and shape of the boundary separating the yielded and unyielded masses of the fluid must be part of the solution of the initial boundary value problem. Several authors have addressed in the past the analysis of the flow of Bingham fluids in conduits and related geometries. The groundbreaking work of Russian researchers in the sixties set the tone for the research direction for decades to come. Safronchik [1–3] and Mosolov and Mjasnikov [4–6] conducted fundamental investigations on the propagation and the location of the yield surface and its properties and the plug and dead regions in the flow, respectively. The channel flow of a Bingham plastic with a given initial velocity distribution and a time dependent pressure gradient is investigated in [1] to determine the subsequent velocity field and the location of the yield surface. A highly complicated equation for the velocity dependent on the location of the interface between the plug zone and the flowing mass of the fluid is Corresponding author at: Centro de Investigación en Creatividad y Educación Superior & Departmento de Ingeniería Mechánica, Universidad de Santiago de Chile, Santiago, Chile. E-mail addresses: [email protected], [email protected] (D.A. Siginer). http://dx.doi.org/10.1016/j.ijnonlinmec.2016.09.012 Received 19 September 2016; Accepted 23 September 2016 Available online 04 October 2016 0020-7462/ © 2016 Elsevier Ltd. All rights reserved. International Journal of Non–Linear Mechanics 88 (2017) 1–10 M.F. Letelier et al. not solved and no examples are given. The existence and uniqueness of the plug region was shown by Mosolov and Mjasnikov [4] using variational methods. In particular they show exactly and analytically that the flow stops if the Bingham 4 number Bi exceeds a critical Bingham number Bic= 2 + π =1.0603178 with no-slip conditions on the walls when the cross-section is a square. This result, that there is a limiting Bingham number for a given crosssectional tube, was extended later to general cross-sectional tubes by Duvaut and Lions [12], and much later was confirmed through extensive numerical computations by Saramito and Roquet [13] for a square cross-section with no-slip conditions as well as with slip on the walls Roquet and Saramito [14]. In the latter case Bic depends on a dimensionless slip parameter S that quantifies the extent of slip. For S =0.6 the critical Bingham number is determined to be 2 Bic= 2 + π =0.5301589. When S =+∞ fluid sticks to the wall and no-slip condition prevails. The existence of dead regions with a concavity always turned towards the inside of the cross-section was proven analytically by Mosolov and Mjasnikov [5]. This is clearly a very difficult proposition to show numerically as a general statement as is obvious from the numerical results presented by Saramito and Roquet [13] and Roquet and Saramito [14]. The first numerical study of yield stress fluids was done by Fortin [15]. The augmented Lagrangian method framework introduced by Fortin and Glowinski [16] and further developed by Glowinski and LeTallec [17] was used by Huilgol and Panizza [18] to study the flow in an annulus and an L-shaped cross-sectional tube. The flow of a Bingham fluid in a square duct, and in particular the geometry of the plug zones as well as the dead zones at the duct corners was explored numerically by Saramito and Roquet [13] and Roquet and Saramito [14]. In the earlier paper [13] the fully developed Poiseuille flow of a yield stress fluid in a square cross-section was studied with no-slip condition at the walls using the augmented Lagrangian method framework coupled to a mixed anisotropic auto-adaptive finite element method. In the later paper [14] the Poiseuille flow of a yield stress fluid in a square cross-section with slip yield boundary condition was considered numerically using the same approach introduced by the authors previously. The consideration of slip is important as it frequently occurs in the flow of two-phase systems such as suspensions, emulsions and in industrial viscoplastic flow problems such as concrete pumping, and it appears to be more pronounced when the material has a yield stress property such as bio-solids and pastes. Steady flow of Bingham fluids in narrow eccentric annuli was investigated both analytically and numerically by Walton and Bittleston [19], and conditions for the existence of plug zones and quasi plug zones were discussed in the context of the flow of pastes and suspensions in complex geometries with industrial applications. Wachs [20] also studied the problem for a wide range of the relevant parameters using numerical methods. The present authors [21,22] studied viscoplastic flows in a variety of non-circular tubes both analytically and numerically relating geometric and flow variables to predict velocity distribution, rate of flow, and energy dissipation. Recent work that contributes insights into the physics of the flow of viscoplastic fluids and the complexities associated with it, as constitutively characterized by the Bingham model, addresses thermal effects, Akram et al. [23] and Turan et al. [24], and the Lattice Boltzmann method applications in complex geometries, Tang et al. [25]. A major difficulty with the flow of a Bingham fluid in complex geometries is the existence of spatial discontinuities in kinematic and dynamic variables at the interface between regions undergoing deformation and the plug and stagnant (or dead) zones. To the best knowledge of the present authors, no published analytical work exists shedding light onto the general behavior of flows of Bingham fluids in non-circular tubes of arbitrary cross-section. The kinematics and dynamics of the steady, developed and isothermal flow of Bingham fluids in longitudinally constant, non-circular cross- derived, and a non-linear integral equation for the position of the interface is presented. The unsteady motion of a Bingham fluid both in bounded and unbounded domains is investigated in [2], specifically the unsteady flow field generated by a rotating cylinder and the flow field in a Couette device, respectively. A non-linear integral equation to determine the location of the yield surface is derived but left unsolved. The unsteady flow of a Bingham fluid in a circular tube is examined in [3] and again an unsolved complicated equation is presented to determine the location of the yield surface. These equations can be solved albeit with considerable difficulty; however whether or not their predictions lead to the correct physical solution, the speed of propagation and position of the yield surface, is an open issue as can be seen from the discussion in the papers of Huilgol [7,8]. For instance the attempt by Huilgol [7] to solve the non-linear integral equation derived by Safronchik [1] to determine the speed and location of the yield surface in the Rayleigh problem for a Bingham fluid led to the nonexistence of a solution. Huilgol states that the reasons behind this nonphysical result are the importance of the homogeneous and nonhomogeneous boundary conditions imposed on the velocity field at the yield surface, and points out that a qualitative understanding of the existence or non-existence of moving yield surfaces in any flow of yield stress fluids requires a deeper insight into the right conditions to be imposed at the interface. Conditions for the existence of the plug zones were investigated by Huilgol [8] who showed in particular the singularity of the yield surface across which the velocity, the acceleration and the velocity gradient as well as the shear stress, its time derivative and its gradient are all continuous, but the time derivative of the acceleration, the spatial gradient of the acceleration, the second gradient of the velocity and the corresponding temporal and spatial gradients of second order of the shear stress all undergo jumps. The pioneering work of Safronchik [1–3] and the work of Huilgol [7,8] shows how difficult it is to determine the location and speed of propagation of a yield surface even in simple geometries. The complexity of the behavior of the yield surface was further brought to light by Glowinski [9] who showed that the yield surface may move laterally with a finite speed in the pressure gradient driven flow of a Bingham fluid, the rigid core in the center gradually becoming larger, the yield surface expanding with a finite speed of propagation, decelerating the fluid and eventually choking off the flow and bringing it to rest. Specifically he proved that the flow of a Bingham fluid in a pipe of arbitrary cross-section comes to a halt in a finite amount of time if the pressure gradient drops suddenly below a critical value needed to overcome the effect of the yield stress in contrast to the behavior of a Newtonian fluid which comes to rest in an infinite amount of time when the pressure gradient is removed suddenly. Sekimoto [10,11] made important contributions in the early nineties to the determination of the propagation speed and the location of the yield surface. He finds in [10] a similarity solution and derives an equation for the location of the yield surface in the case of an existing steady simple shear flow in the semi-infinite region over a flat plate which becomes unsteady due to a sudden reduction in the shear stress on the boundary to a value below the yield stress of the Bingham fluid. He shows that the yield surface propagates from the flat plate boundary into the fluid and derives an evolution equation for the location of the yield surface at subsequent times; however he does not solve the equation either. The lateral movement of a yield surface in a shearing flow is considered in [11]. He correctly assumes that the lateral motion of the yield surface obeys the diffusion equation, ρ ∂u ∂τxy = ∂t ∂y and asserts that the gradient of the shear stress τxy is continuous across the yield surface provided certain continuity conditions are met the continuity of the local acceleration among them. He presents an evolution equation for the location of the yield surface; however it is 2 International Journal of Non–Linear Mechanics 88 (2017) 1–10 M.F. Letelier et al. terms of stress is written as, sectional tubes are investigated in this paper. An analytical method previously developed by the authors [21,26–28] for modeling flow in tubes of arbitrary contours is extended to viscoplastic flow. The method is based upon a one-to-one mapping of a reference circle representing a tube of circular cross-section into cross-sections of arbitrary shape. This approach makes it possible to solve the equations of motion in a great variety of cross-sectional shapes. Tube contours can be selected from a continuous spectrum of shapes. In particular, the method is applied to the analysis of central plug zones, stagnant zones at sharp corners, shear stress distribution and energy dissipation, and rate of flow variation with yield stress. Results presented are partially consistent with previously published investigations. Certain analytical anomalies that deserve close attention are described with some detail, since they point out aspects of the physics of Bingham flows that may have been overlooked so far. Although Bingham constitutive equation is widely used to model the flow behavior of viscoplastic fluids to predict the evolution of the flow in many areas of application it may fall short of expectations under certain prevailing conditions as it includes only two material parameters, namely the yield stress and viscosity, which may restrict the use of that model when either the constitution of the fluid is more complex or the flow characteristics demand more involved predictive capacity. In these cases the Herschel-Bulkley (H-B) constitutive equation, an extension of the Bingham model with a third material parameter, i.e. a power index, can be used. The range of values of the power index of the H-B model greater and lesser than unity determines the extent of the shear-thinning and shear-thickening behavior, respectively. H-B constitutive model has been extensively applied successfully, in particular recently, to a variety of problems such as the dispersion of a solute in an annular tube for modeling blood flow in an artery affected by a catheter, Ramana & Sarojamma [29], non-Newtonian flows with a stratified interface, Moyers-Gonzalez et al. [30], and peristaltic pumping in a channel, Vajravelu et al. [31] among others. The analytical method introduced in this paper to predict the flow field of a Bingham plastic is used to determine the flow field of an H-B fluid in tubes of arbitrary cross-sectional geometries. The analysis is valid for a wide spectrum of non-circular cross sectional shapes. As a specific example the influence of shear-thinning and shear-thickening on the plug zones and stagnant zones has been investigated in the flow field in a tube with equilateral triangular cross-section and a comparative study with the predictions of the Bingham model is presented. The paper is organized as follows: The flow of Bingham plastics and Hershel-Bulkley fluids in straight tubes of arbitrary cross-section is analyzed in Sections 2 and 3, respectively. Results are discussed in Section 4 and Conclusions are presented in Section 5. ∂τrz τrz 1 ∂τθz ∂P + + =− ∂r r r ∂θ ∂z in which P is the pressure. The analysis of the flow structure can be enhanced by the use of natural coordinates, Letelier and Siginer [21]. If n is a curvilinear coordinate normal to the flow isovels (equal velocity lines), then (1, 4) transform into, τnz=Bi − dw , τnz≥Bi dn (5) dτnz τnz ∂P + =− ∂z dn ρ (6) in which τnz is the normal axial shear stress and ρ is the radius of curvature of isovels. Flow analysis in a range of tube cross-sectional shapes can be achieved by imposing a boundary condition through a shape factor G defined as, G=1 − r 2+εr αsinαθ (7) where (ε, α ) are mapping parameters of a one-to-one and continuous mapping to map the base circumference into a continuous spectrum of tube contours. ε is chosen as the perturbation parameter in the analysis presented in this paper. The parameter α must be an integer number in order to get closed shapes for the mapped tube contour. Combinations of ε and α yield many different contour shapes forG = 0 while simultaneously satisfying the no-slip boundary condition as it will become apparent in the analysis below. This shape factor can be made more general by adding more terms to the mapping (7). The axial velocity w is expressed as an asymptotic expansion, w=G (r , θ )[ f0+εf1+ε 2f2 +O (ε 3)] (α −2)/2 2 ⎛ α −2 ⎞ εc= ⎜ ⎟ α⎝ α ⎠ (9) Introducing the expansion of the velocity field w into asymptotic series, 2. Analysis: bingham fluids and substituting (10) into (1) we obtain, For parallel flow in tubes, Bingham constitutive equation, in cylindrical coordinates, is given by, τrz=Bi − (1) τrz≥Bi , τθz≥Bi (2) ⎛ ∂w1 ⎞2 ⎤ 2 Bi ∂wo ⎡ ∂w1 ⎤ ⎡ ∂w2 −⎢ + 2 ⎟ ⎥ ε + O (ε 3 ) ⎜ ⎥ ε−⎢ ∂r ⎣ ∂r ⎦ ⎢⎣ ∂r 2r (∂wo /∂r )2 ⎝ ∂θ ⎠ ⎥⎦ ⎞ 1 ⎛ ∂w τθz=− ⎜ o −N ⎟ ⎠ r ⎝ ∂θ ∂w1 ∂θ ∂wo ∂θ (10) (11) ⎡ ⎛ ∂w1 ⎞ ⎛ ∂w1 ⎞ N 2 ∂w2 ⎤ 2 ⎥ ε + O (ε 3 ) ε−⎢ ⎟− ⎟⎜ ⎜ ⎣ r (∂wo /∂θ )2 ⎝ ∂θ ⎠ ⎝ ∂r ⎠ (∂wo /∂θ ) ∂θ ⎦ (12) 1 2 ∂w ⎞ ⎤ 2 ⎡ ⎛ ∂w ⎞2 ⎛ 1 I =⎢ ⎜ ⎟ + ⎜ ⎟⎥ ⎢⎣ ⎝ ∂r ⎠ ⎝ r ∂θ ⎠ ⎥⎦ (8) in which the functions fi are unknowns to be determined by substituting Eq. (8) into Eq. (4). The structure of (8) ensures that asymptotic conditions for ε = 0 and for Bi = 0 are met. It should be noted that the mapping parameter (also the perturbation parameter) ε has maximum allowable values εc for a given value of α . If ε > εc the mapped contour is no longer closed. The limiting value εc of ε is given by the following expression when α > 2 , Siginer and Letelier [26,27], w0=w0 +εw1+ε 2w2+O (ε 3) ⎛ ⎛ Bi ⎞ ∂w Bi ⎞ 1 ∂w , τθz=−⎜1 + τrz=−⎜1 + ⎟ ⎟ ⎝ ⎝ I ⎠ ∂r I ⎠ r ∂θ (4) The zero order velocity is the well-known Bingham formula for flow in round tubes, i.e. (3) w0=1 − r 2−Bi (1 − r ) where τrz and τrθ are the axial shear stress components, Bi is the Bingham number a dimensionless yield stress, and I defined in (3) is the second invariant of the rate of deformation tensor. Scale factors a (a reference radius) for the radial coordinate r , w0 for the axial velocity w, ηw0 / a for all shear stress components, and ηw0 / a2 for the axial pressure gradient, in which η is the fluid viscosity have been introduced to render all variables dimensionless. The momentum equation in (13) Thus f0 =1 − Bi 1+r (14) The dimensionless pressure gradient is (−4). Substitution of (11– 14) in (4) gives, 3 International Journal of Non–Linear Mechanics 88 (2017) 1–10 M.F. Letelier et al. ⎛ Bi ⎞ ⎛ ∂ 2w 1 ∂w1 ⎞ ∂ 2w1 r ⎜r − ⎟ ⎜ 21 + =0 ⎟+ ⎝ 2 ⎠ ⎝ ∂r r ∂r ⎠ ∂θ 2 where K (r , θ ) is a complex function that involves terms in w0 and w1. Eq. (21) can be solved analytically and hence f2 can be found in a similar fashion to f0 and f1. Computing f2 does not change results qualitatively as compared to the case where only f0 and f1 are considered. Analytical expressions for fi ,i = 3, …, M can be determined in a similar way but make only increasingly small quantitative contributions to the values of the flow variables, and do not change at all the qualitative outcomes as would have been implied by the asymptotic nature of the expansions. Energy dissipation can be computed through the dissipation function ϕ . In natural coordinates ϕ is given by (15) The mathematical structure of the longitudinal velocity field w implies that the tube boundary can be specified in an infinite number of different forms. The range of the spectrum of mapped closed tube contours may be substantially extended if an infinite number of additional mapping terms of the form εi r αi sin αi θ , εi r αi cos αi θ are incorporated into the shape factor G defined in (7). That makes it possible to describe very complex shapes depending on the number of additional mapping terms included. Shape factor G describes exact ellipses for α = 2 and any ε , an exact equilateral triangle for the limiting case of ε = εc=0.3849 , and regular figures of number of sides equal to α in which sides are slightly curved when ε = εc as in (9) if only one mapping term is considered in (7). In the present analysis two families of shapes are considered; shapes that evolve from a circle to an equilateral triangle (α = 3,0 ≤ ε ≤ εc=0.3849) and shapes that evolve from a circle to a square (α = 4,0 ≤ ε ≤ εc=0.25) with slightly curved sides when ε = εc . These tube contours were selected because they are the contours with the largest departure from the circle the base contour, and thus, are more likely to bring out the dynamic characteristics of Bingham fluid flows. To solve (15) we set, w1=Rα (r ) sinαθ ϕ = τn Clearly ϕ is a function only of the normal shear stress at any point with the yield stress, the Bingham number as a parameter, that is, ϕ=ϕ (τn; Bi ) 3. Analysis: Hershel-Bulkley (H-B) fluids The constitutive equation of an H-B fluid for parallel, 2-D flow, in dimensionless variables, reduces to, (16) ⎛ ⎛ Bi ⎞ 1 ∂w Bi ⎞ ∂w , τθz=−⎜I n −1+ ⎟ τrz=−⎜I n −1+ ⎟ ⎝ ⎝ I ⎠ r ∂θ I ⎠ ∂r where α is the mapping parameter in (8). We make a change of 2r dependent variable r= Bi and transform (15) into r (1 − r ) d 2Rα (r ) dR (r ) 2 +(1 − r ) α +α Rα (r )=0 dr 2 dr The first hypergeometric function F12 (a1, a2; b1; r ) corresponding to p=2, q =1 is the solution of (17). It is also known as “the” hypergeometric equation or Gauss’s hypergeometric function, and it arises frequently in physical problems. The solution of (17) reads, ⎛ 2⎞ 2 − Bi −1 ⎟ 2 F1 (a, b ; c; r ) 2 F1 ⎜a , b ; c ; ⎝ 2 Bi ⎠ (a , b ; c ; r ) , a = − α , b = α , w0=K [(Y (1))−Y (r )]; K = c=1 Y (r )=(2r − Bi ) β ; β = 2 F1 (− α , α ; 1; r ) is reduced to polynomials for positive, integer values of α . Evaluating the solution for α=3 and α=4 corresponding to the equilateral triangular and the square cross-sections we obtain, respectively, R3 and R 4 , Bi 3 9 2 9 R3 = − + Bi r − Bi r 2 + r 3 80 40 10 f0 = (19) (25) (26) ∂ 2w ∂w ∂ 2w1 +[2(2n −1) r −nBi ] 1 +2 21 =0 ∂θ ∂r ∂r 2 (27) The function f1 in the asymptotic expansion (8) is computed through the solution of Eq. (27). The solution of Eq. (27) for the description of the H-B fluid flow in a tube with arbitrary cross-section is sought by prescribing And for the equilateral triangular cross-section α=3 we determine, w1 (r , θ )=Rα (r )sin αθ (20) (28) where Rα (r ) is a function of the radial coordinate r, and α is the parameter appearing in the shape factor Eq. (1). Substitution of Eq. (28) into Eq. (27) yields the following equation for Rα (r ) The second order velocity w2 can be computed from the differential equation of second order derived from (4), ⎛ Bi ⎞ ⎛ ∂ 2w 1 ∂w2 ⎞ ∂ 2w2 r ⎜r − ⎟ ⎜ 22 + =K (r , θ ) ⎟+ ⎝ 2 ⎠ ⎝ ∂r r ∂r ⎠ ∂θ 2 w0 1−r 2 nr (2r −Bi ) Rα (r ) − r αf0 sinαθ 1−r 2 (2 − Bi )(1 + r ) R3 (r ) − 2 R3 (1)(1 + r − Bi ) r 3 sin3θ 2R3 (1)(1 + r )(1 − r 2 ) 1+n n (24) Shear stress components can be linearized through a combination of Eqs. (22–25) and substituted into the linear momentum balance Eq. (4) following a similar procedure as in Section 2. The resulting equation for w1 in the asymptotic expansion (10) is given by, (18) Equating terms of O(ε) in (8) and (10), and imposing the regularity condition that f1 be continuous at r = 1 we find, f1 = n 2(1+n ) The comparison of O(1) terms in Eq. (8) and Eq. (10) yields, R 4 = 0. 00179Bi 4 − 0. 05714Bi 3r + 0. 42857Bi 2r 2 − 1. 14286Bi r 3 + r 4 f1 = (23) τrz andτθz are the components of the extra-stress tensor, and Bi and n represent the dimensionless yield stress, the Bingham number, and the power index, respectively. The range of the power index specifies the shear-thinning and shear-thickening capability of the H-B fluid. The axial velocity field w(r , θ ) is expanded in an asymptotic series similar to the expansion (10) in Section 2 in terms of the small one-to-one mapping parameter ε. In this case w0 in the asymptotic expansion (10) is the well-known velocity field of an H-B fluid in a circular tube, i.e. Fqp (a1, ….,ap; b1, ….,bq; r ) 2 F1 (22) 1 ⎡ ⎛ ∂w ⎞2 ⎛ 1 ∂w ⎞2 ⎤ 2 I =⎢ ⎜ ⎟ + ⎜ ⎟⎥ ⎢⎣ ⎝ ∂r ⎠ ⎝ r ∂θ ⎠ ⎦⎥ (17) This is the hypergeometric equation of the second kind the solution of which is given in terms of the generalized hypergeometric function R α (r ) = dw =τn (τn−Bi ) dn nr (2r −Bi ) (21) 4 d 2Rα (r ) dR (r ) +[2(2n −1) r −nBi ] α −2α 2Rα (r )=0 dr 2 dr (29) International Journal of Non–Linear Mechanics 88 (2017) 1–10 M.F. Letelier et al. Making the same type of dependent variable change r =2r / Bi as in Section 2 reduces (29) to the following Hypergeometric equation, r (1 − r ) ⎛ 2n −1 ⎞ ⎤ dRα (r ) α 2 d 2Rα (r ) ⎡ + ⎢1 − ⎜ + Rα (r )=0 ⎟r ⎥ ⎝ n ⎠ ⎦ dr dr 2 n ⎣ (30) The regular solution to Eq. (13) is given in terms of the Hypergeometric function of the first kind F12 (a, b;c;r ) Rα (r ) = γF12 (a, b; c; r ) a=− 1−n +∆ 1−n −∆ ; b=− ; c=1;∆ = 2n 2n (1−n )2 +4α 2n where γ is a constant whose value is determined by imposing a regularity condition at r = 1, γ= (2−Bi )1/ n 2 α = 3; ε = 0.3849, Bi=0.3 Fig. 2. Normal shear stress distribution for τn=0.3,0.5,0.7,0.9,1.1,1.3,1.5,1.7,1.9,2.1and 2.3 from center to wall. 2 F12 (a, b; c; Bi ) It is worthwhile to remark that for Bi = 0 the function Rα (r ) reduces to the O(ε) contribution for the purely power-law fluid solution, i.e. Bi →0 Rα ⎯⎯⎯⎯⎯→ 2 1− n −n n r Similarly, for n=0, Rα (r ) collapses onto the O(ε) solution for the purely Bingham flow as previously found for the triangular and square cross-sections, namely, for α=3, n →0 ⎛ Bi3 9Bi2 9Bi 2 3⎞ R3 ⎯⎯⎯⎯→ ⎜ − + r− r +r ⎟ ⎝ 80 ⎠ 40 10 Finally, comparison of O(ε) terms in Eq. (8) and Eq. (10), and using Eq. (28) yields f1 as f1 = Rα (r )−r αf0 sin αθ 1−r 2 Higher order terms can be determined in like manner. Fig. 3. Equal velocity lines for α = 3; ε = 0.3849, Bi=1.65, w=0, 0.007, 0.014 and 0.02 from border to center. Plug zone is in light gray color. Stagnant zones are in dark gray color in the corners. 4. Results Setting G(r, θ) =0 and varying the values of the mapping parameters ε and α yields a great variety of cross-sectional contours. The values ε = 0.3849 and α = 3 leads to an exact equilateral triangle; the geometry corresponding to ε = 0.22 and α = 4 is a square with rounded corners. The velocity field and normal shear stress distribution for α = 3 are shown in Figs. 1–9 and Figs. 10–11 show the velocity field for α = 4 . Figs. 1–9 show the evolution of the plug and stagnant zones for the equilateral triangular cross-section for three values of the yield stress. The boundary of the plug zones is defined as the maximum velocity equal velocity line, and the stagnant regions at the corners are limited by the walls and the line on which w = 0 . This definition for the Fig. 4. Normal shear stress distribution for α = 3; ε = 0.3849, Bi=1.65 τn=1.683,1.8,1.9 and 2.0 from center to wall. plug boundary makes physical sense since a solid region must move with constant velocity. Fig. 7 illustrates a 3-D plot of the mathematical predictions of (8), where negative velocity contours at the center and inside the stagnant regions have been eliminated. This full mathematical configuration is complex and is only partly realistic. Predicted local maxima, as shown in Fig. 8 for Bi = 1.7, cannot exist because they violate the condition stated in (2), that is, shear cannot be smaller than the yield stress, which would be the case if there are zones of positive velocity gradient in the flow. Therefore the only realistic interpretation is that the Bingham model is valid from the tube wall up to the Fig. 1. Equal velocity lines for α = 3; ε = 0.3849, Bi = 0.3, w = 0,0.15,0.3,0.45,0.6and0.722 from border to center. Plug zone is in gray color. 5 International Journal of Non–Linear Mechanics 88 (2017) 1–10 M.F. Letelier et al. Fig. 5. Equal velocity lines for α = 3; ε = 0.3849, Bi =1.7 w=0, 0.0111 from border to center. Plug zone is in light gray color. Stagnant zones are in dark gray color in the corners. Fig. 8. Equal velocity lines corresponding to Fig. 7. Fig. 9. Plug and stagnant zones for α = 3, ε = 0.375, Bi = 1.7, w = 0.0120 . Plug zone is in light gray color. Stagnant zones are in dark gray color. Fig. 6. Normal shear stress distribution for α = 3; ε = 0.3849, Bi =1.7; τn=1.736,1.8,1.9 and 2.0 from center to wall. Fig. 7. 3-D plot of velocity field for α = 3; ε = 0.3849, Bi=1.7. continuous curve of maximum velocity. However, in these cases there would be a velocity gradient discontinuity in parts of the plug zone contour. Physically this means, as can be appreciated from Figs. 2, 4 and 6, that the normal shear stress is not constant in the plug boundary, as it is for smaller . Thus the discontinuity is of mathematical nature but not physical. The reliability of these analytical results have been confirmed through the following validation procedure; second order terms in ε have been considered to show that they have little effect on the results, and the calculations have been extended to contours with much smaller value for ε than for α = 3 and ε = 0.3849 : α = 4 (ε = 0.204), α = 5 (εc=0.186) Fig. 10. Plug and stagnant zones for α = 4, ε = 0.25, Bi = 1.8, w = 0.00297 . Plug zone is in light gray color. Stagnant zones are in dark gray color. and α = 6(εc=0.148). In all these cases same type of flow configurations prevails. The work by Saramito and Roquet [13,14] shows concave stagnant zones at the corner, whereas the present analysis predicts convex ones, using the same definition criteria. The present method provides a mathematical argument that backs the convex shape prediction. Eq. (6) shows that at the corner apex τn is zero since the radius of curvature is zero. Deformation starts when τn≥Bi . Figs. 2, 4 and 6 show that close to 6 International Journal of Non–Linear Mechanics 88 (2017) 1–10 M.F. Letelier et al. computed for α = 3 and ε = 0.3849 and three values of the dimensionless yield stress the Bingham number Bi, Bi 0.5 1.0 1.5 Q 1, 171 0, 620 0, 170 A 100% increase in the dimensionless yield stress Bi from Bi = 0.5 to Bi = 1.0 decreases Q by 47%. In all these computations the pressure gradient has been kept fixed at a value of −4. Other values of the pressure gradient do not change the general outcomes of this analysis. In steady flow axial pressure forces are balanced out by axial shear stress for any fluid element. In circular tubes constant normal shear stress curves are also circular, so that all circular concentric elements develop a constant shear stress on their limiting surfaces. This is not the case for non-circular tubes, except when the dimensionless yield stress Bi and the radius of the fluid element are small, see Fig. 1. Figs. 2, 4 and 6 depict the curves of constant normal shear stress which are directly related to the dissipation function through Eq. (2). The maximum normal shear stress occurs at the center of the walls of the triangle, and the minimum τn=0 occurs at the apex of the corners and at the tube center. Greater normal shear stress means greater normal velocity gradients implying that velocity increases faster closer to the middle of the walls of the triangle. Knowledge of energy dissipation distribution inside the flow region can help improve conduit design and help predict head loss in flows in complex geometries. Next the general expression obtained for the velocity field of an H-B fluid is evaluated for the equilateral triangle case parametrized by the values of α and ε . The results for the square with rounded corners are not much different qualitatively with respect to the presented case and are not displayed. Equal velocity curves (isovels) and velocity profiles are shown for the triangular cross-section corresponding to α = 3 and ε = 0.3849 . The yield stress takes in each figure the following values: Bi∈{0.0,0.5,1.0,1.3}. Within each figure, the behavior of a shear-thin- Fig. 11. Plug and stagnant zones for α = 4, ε = 0.204, Bi = 1.8, w = 0.00576 . Plug zone is in light gray color. Stagnant zones are in dark gray color. the corners iso-τn curves are convex, a result that can be very accurately computed for smaller Bi and is exact for Bi = 0. The main argument is that the normal shear stress curves at the corners are convex, starting with Newtonian fluids, and at the border of the dead zones the normal shear stress should be equal to the yield stress, and that is what the analysis shows. The effect of sharp corners on the flow is illustrated in Figs. 9–11, in which sharp edges have been smoothed by taking a smaller value of ε with everything else equal. In these cases the stagnant zones decrease significantly as the plug velocity increases. Rate of flow Q = ∫ wrdrdθ s can be computed once the geometry of isovels and solid zones is known. S is a piece-wise area composed of a central constant velocity portion and a complex-shaped annulus in which the velocity varies. Stagnant zones are not considered. Three values of Q have been Fig. 12. Isovels and velocity profiles along a symmetry axis for a fixed value of the Bingham number Bi=0 and various values of the power-law index: (a) and (d): n = 0.5 ; (b) and (e): n = 1.0 ; (c) and (f): n = 1.5. 7 International Journal of Non–Linear Mechanics 88 (2017) 1–10 M.F. Letelier et al. Fig. 13. Isovels and velocity profile along a symmetry axis for a fixed value of the Bingham number Bi=0.5 and various values of the power-law index: (a) and (d): n = 0.5 ; (b) and (e): n = 1.0 ; (c) and (f): n = 1.5. ning fluid (n = 0.5), purely plastic fluid (n = 1), and shear-thickening fluid (n = 1.5) is presented in the 1st, 2nd and 3rd columns, respectively. The corresponding velocity profile is drawn along a symmetry axis of the cross-section. For Bi ≠ 0 , plug-zones (PZ) are identified by the contour given by the maximum velocity isovel. Plug-zones at the cross-sectional center and stagnant zones (SZ) at the borders are shown (if any) as shaded areas with corresponding velocities shown in the velocity profiles. Fig. 12 depicts the behavior of a purely power-index fluid, Bi=0 , where neither SZ or PZ are present. When the yield stress value increases up to Bi=0.5 as in Fig. 13 the situation changes and SZ and PZ appear as a function of n. Fig. 13 shows that for n=0.5 there exist SZ but no PZ, while for n=1 and n=1.5 the opposite situation is found. Figs. 14 and 15 show that this difference becomes more pronounced as Bi takes higher values. In each figure, instead of a PZ for n=0.5 the central zone of the cross-section holds comparatively high velocities, with a cusp at r = 0 which gets sharper as Bi increases. Hence, even if the fluid becomes stagnant near the edges of the cross-section it is rapidly flowing at the center. Conversely, when n=1.5 the opposite occurs, i.e. PZ are formed but no SZ appears for the values of the yield stress analyzed. Moreover, for any fixed value of n , as the Bingham number Bi increases, the velocity and hence, rate of flow decreases in line with plastic behavior. cross-sectional geometry is determined by the shape factor method, which maps the circular base contour onto a family of shapes characterized by two parameters, α and ε . The exact equilateral triangular cross-section (α = 3, ε = 0.3849 ) has been chosen to illustrate explicit results. However it is relevant to point out that the velocity field can be evaluated for any arbitrary cross-sectional contour via the shape factor method. Viscoplastic flows in non-circular tubes develop zones of undeformed fluid called plug zones and if the tube contour has sharp corners, stagnant zones appear in the vicinity of the apex of the corner. We show that the Bingham model may mathematically predict velocity fields that are in part not realistic and need physical interpretation. This is especially so when the plug zones have non-circular boundaries. In this case the longitudinal shear stress is not constant at the boundary of the plug zone, and it is equal or greater than the yield stress. The analytical method presented yields fast results for a wide variety of contours, thus allowing the comparison and quantification of the effect of different contours on the kinematic and dynamic variables of the flow. We show in particular that a radius of curvature at the corners no matter how small greatly diminishes the stagnant zones and increases the rate of flow, which means that even a very small departure from sharpness at the corners produces this result. In addition we also predict that stagnant areas at the corners of the cross-section are convex in contrast to the findings of previously published work, Saramito and Roquet [13] and Roquet and Saramito [14], and we provide a rigorous mathematical argument in support of our predictions. The rate of flow and energy dissipation distribution in the tube cross-section has been explored. Energy dissipation, an important variable in the design of conduits for viscoplastic materials, is closely related to normal shear stress distribution, which can be very accurately determined via the analytical approach presented in this paper. Flow of a Herschel-Bulkley fluid in tubes of non-circular crosssection has been investigated analytically using the same approach as 5. Conclusions The extension of a versatile analytical method previously used by the authors for steady flow in non-circular tubes of viscoelastic fluids to the flow of Bingham fluids has provided several novel insights concerning the flow dynamics. This method was used to study the secondary flows and heat transfer in flows of viscoelastic fluids leading to novel analytical predictions, Siginer & Letelier [26,27]. In this paper the analysis is focused mainly on an equilateral triangular tube crosssection, and also on the square cross-sectional tube to some extent. The 8 International Journal of Non–Linear Mechanics 88 (2017) 1–10 M.F. Letelier et al. Fig. 14. Isovels and velocity profile along a symmetry axis for a fixed value of the Bingham number Bi=1.0 and various values of the power-law index: (a) and (d): n = 0.5 . (b) and (e): n = 1.0 ; (c) and (f): n = 1.5. for Bingham fluids. In general the results are in agreement with the standard behavior of power-index fluids in circular pipes. For any value of the Bingham number Bi and power index n > 1 that is when the fluid is shear-thickening, velocity increases when compared to the case of the purely plastic Bingham fluid (n = 1). For shear-thinning fluids, that is when n < 1, the opposite behavior occurs, that is overall the velocity Fig. 15. Isovels and velocity profile along a symmetry axis for a fixed value of the Bingham number Bi=1.3 and various values of the power-law index: (a) and (d): n = 0.5 ; (b) and (e): n = 1.0 ; (c) and (f): n = 1.5. 9 International Journal of Non–Linear Mechanics 88 (2017) 1–10 M.F. Letelier et al. 1972. [13] P. Saramito, N. Roquet, An adaptive finite element method for viscoplastic fluid flows in pipes, Comput. Methods Appl. Mech. Eng. 190 (2001) 5391–5412. [14] N. Roquet, P. Saramito, An adaptive finite element method for viscoplastic flows in a square pipe with stick-slip at the wall, J. Non-New. Fluid Mech. 155 (2008) 101–115. [15] M., Fortin, Calcul numérique des écoulements des fluids de Bingham et des fluides Newtoniens incompressibles par la méthode des éléments finis, thèse, Paris VI, 1972. [16] M. Fortin, R. Glowinski, Méthodes de Lagrangian augmenté. Applications à la résolution numérique de problèmes aux limites, Méthodes Mathématiques de l′Informatique, Dunod, 1982. [17] R., Glowinski and P., Le-Tallec, Augmented Lagrangian and operator-splitting methods in non-linear mechanics, Stud. Appl. Math. Soc. Ind. Appl. Math. 1989. [18] R.R. Huilgol, M.P. Panizza, On the determination of the plug flow region in Bingham fluids through the application of variational inequalities, J. Non-New. Fluid Mech. 58 (1995) 207–217. [19] I. Walton, S. Bittleston, The axial flow of a Bingham plastic in a narrow eccentric annulus, J. Fluid Mech. 222 (1991) 39–60. [20] A. Wachs, Numerical simulation of steady Bingham flow through an eccentric annular cross-section by distributed Lagrange multiplier fictitious domain and augmented Lagrangian methods, J. Non-New. Fluid Mech. 142 (2007) 183–198. [21] M.F. Letelier, D.A. Siginer, On the flow of a class of viscoinelastic-viscoplastic fluids in tubes of non-circular contour, Int. J. Eng. Sci. 45 (2007) 873–881. [22] M. F. Letelier, F. Godoy, and C. Rosas, On the existence conditions for plug zones in plastic flow in tubes of non-circular cross-section, in: Proceedings of the 14th PanAmerican Congress of Applied Mechanics, Santiago, Chile, 2014. [23] S. Akram, S. Nadeem, A. Hussain, Effects of heat and mass transfer on peristaltic flow of a Bingham fluid in the presence of inclined magnetic field and channel with different wave forms, J. Magn. Magn. Mater. 362 (2014) 184–192. [24] O. Turan, N. Chakraborty, R. Poole, Laminar natural convection of Bingham fluids in a square enclosure with differentially heated side walls, J. Non-New. Fluid Mech. 154 (2010) 901–913. [25] G. Tang, S. Wang, P. Ye, W. Tao, Bingham fluid simulation with the incompressible lattice Boltzmann model, J. Non-New. Fluid Mech. 166 (2011) 145–151. [26] D.A. Siginer, M.F. Letelier, Laminar flow of non-linear viscoelastic fluids in straight tubes of arbitrary contour, Int. J. Heat Mass Transf. 54 (2011) 2188–2202. [27] D.A. Siginer, M.F. Letelier, Heat transfer asymptote in laminar flow of non-linear viscoelastic fluids in straight non-circular tubes, Int. J. Eng. Sci. 48 (2010) 1544–1562. [28] M.F. Letelier, J. Stockle, A shape-factor method for modelling parallel and axiallyvarying flow in tubes and channels of complex cross-section shapes, Biomed. Sci. Eng. Technol. INTECH Croat. (2011) 469–486. [29] B. Ramana, G. Sarojamma, Exact analysis of unsteady convective diffusion in a Herschel-Bulkley fluid in an annular pipe, Int. J. Math. Comput. 3 (1) (2013) 15–26. [30] M. Moyers-Gonzalez, K. Alba, S.M. Taghavi, I.A. Frigaard, A semi-analytical closure approximation for pipe flows of two Herschel-Bulkley fluids with a stratified interface, J. Non-New. Mech. 193 (2013) 49–67. [31] K. Vajravelu, S. Screenadh, V., Ramesh Babu, Peristaltic pumping of a HerschelBulkley fluid in a channel, Appl. Math. Comp. 169 (1) (2005) 726–735. decreases when compared to the purely plastic Bingham fluid, but locally increases in the central region. When the velocity is normalized with the average velocity, the opposite trend is true. A particularly interesting result is that the plug-zone size and geometry are functions of the power law index n. For a fixed value of the Bingham number Bi the plug-zone shape departs from a circle as n moves away from unity, and can become completely absent for shear thinning fluids. The stagnant zones also depend on n . In the present analysis they are present only in the velocity distributions for n = 0.5, but are absent in the other two cases analyzed. Acknowledgements This work was supported by FONDECYT-Chile through grant number 1130346, and by DICYT at the Universidad de Santiago de Chile. References [1] A.I. Safronchik, Non-steady flows of a visco-plastic material between parallel walls, PMM, 23, 925–935, J. Appl. Math. Mech. 23 (1959) 1314–1327. [2] A.I. Safronchik, Rotation of a cylinder with a variable angular velocity in a viscoplastic medium, PMM, 23, 1051–1056, J. Appl. Math. Mech. 23 (1959) 1504–1511. [3] A.I. Safronchik, Unsteady flow of a visco-plastic material in a circular tube, PMM, 24, 149–153, J. Appl. Math. Mech. 24 (1960) 200–207. [4] P.P. Mosolov, V.P. Mjasnikov, Variational methods in the theory of viscous-plastic medium, PMM, J. Appl. Math. Mech. 29 (1965) 468–492. [5] P.P. Mosolov, V.P. Mjasnikov, On stagnant flow regions of a visco-plastic medium in pipes, PMM, J. Appl. Math. Mech. 30 (1966) 705–719. [6] P.P. Mosolov, V.P. Mjasnikov, On qualitative singularities of the flow of a viscoplastic medium in pipes, PMM, J. Appl. Math. Mech. 31 (1967) 581–585. [7] R.R. Huilgol, On kinetic conditions affecting the existence and non-existence of a moving yield surface in unsteady unidirectional flows of Bingham fluids, J. NonNew Fluid Mech. 123 (2004) 215–221. [8] R. Huilgol, On the description of the motion of the yield surface in unsteady shearing flows of a Bingham fluid as a jerk wave, J. Non-New. Fluid Mech. 165 (2010) 65–69. [9] R. Glowinski, Sur l′écoulement d′un fluide de Bingham dans une conduite cylindrique, J. Mécanique 13 (1974) 601–621. [10] K. Sekimoto, An exact non-stationary solution of simple shear flow in a Bingham fluid, J. Non-New. Fluid Mech. 39 (1991) 107–113. [11] K. Sekimoto, Motion of the yield surface in a Bingham fluid with a simple-shear geometry, J. Non-New. Fluid Mech. 46 (1993) 219–227. [12] G. Duvaut, J.L. Lions, Les Inéquations en méchanique et en physique, Dunod, 10