On the Flow of Viscoplastic Fluids in Non-Circular Tubes

International Journal of Non-Linear Mechanics 153 (2023) 104408 Contents lists available at ScienceDirect International Journal of Non-Linear Mechanics journal homepage: www.elsevier.com/locate/nlm On the flow of viscoplastic fluids in non-circular tubes Mario F. Letelier, Dennis A. Siginer ∗, Edgardo González Department of Mechanical Engineering, Universidad de Santiago de Chile, Santiago, Chile ARTICLE INFO Keywords: Viscoplastic fluid flow non-circular tube contour Herschel–Bulkley fluid Bingham fluid Plug and stagnation zones ABSTRACT Steady flow of the viscoplastic Bingham and Herschel–Bulkley (H–B) fluids in tubes of non-circular cross-section is investigated analytically. The solution methodology is general in scope, does not put any restrictions on the Bingham number and thus allows the mapping of the flow field for high Bingham numbers in straight tubes with non-circular axially-symmetric but otherwise arbitrary cross-sectional contours. The circular tube contour is mapped onto an arbitrary non-circular contour on which the no-slip condition is satisfied via a one-to-one and continuous mapping. Governing equations are solved for the full spectrum of axially symmetrical crosssectional shapes and a specific example is developed for the viscoplastic H–B fluid flow in a tube with an equilateral triangular cross-section. The shape and the extent of the plug zones centered on the tube axis and the stagnant zones in the corners are predicted for both Bingham and H–B fluids. The effect of the shear rate dependent viscosity of the H–B fluids, leading to either shear-thickening or shear-thinning behavior, on the formation of the plug and stagnation zones is examined. 1. Introduction Modeling of industrial and biological flows often must deal with viscoplastic fluids that appear in contexts such as mining, food, cosmetic, pharmaceutical, paint, coating and construction industries, blood and other biological flows, and in some flows in nature such as mud and lava displacements. Other important applications concern removal of the air bubbles from yield stress fluids cement and concrete to increase the material strength and the need to remove flammable gas bubbles from another yield stress fluid, radioactive sludge for safety purposes. Bingham constitutive equation [1,2] is widely used to model the flow behavior of viscoplastic fluids for its ability to predict qualitatively the evolution of the flow in most of the above-mentioned areas of application. However, Bingham’s model 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. The Herschel–Bulkley constitutive equation [3], an extension of the Bingham model with a third material parameter, i.e. a power index, allows the mathematical modeling of the flow of viscoplastic fluids with shear rate dependent viscosity. The H–B model has been extensively applied to a variety of problems such as the dispersion of a solute in an annular tube for modeling blood flow in an artery, solute dispersion, mass transfer across tube boundary and peristaltic pumping in a channel just to name a few among many others. Perhaps the main difficulty in computations with yield stress fluids is the location of the yield surface that is generally unknown. The jump in the kinematic variables at the yield surface poses a serious problem in numerical computations although for analytical studies it is not so much of an obstacle. A common approach in numerical studies is to adopt a regularization scheme to remove the discontinuity at the yield surface. The regularization methods of Bercovier and Engelman [4] and Papanastasiou [5] are used widely. An extensive background of the research into the flow of viscoplastic fluids is given in Letelier et al. [6,7]. A recent overview of the pending issues with viscoplastic fluids, in particular from the computational point of view among others, can be found in Denn and Bonn [8]. In this paper we investigate analytically the flow field of the H–B fluids in tubes of arbitrary cross-sectional geometry and present as a specific example of the general solution methodology a comparative study of the flow fields of the purely shear thinning and thickening power law fluids with purely viscoplastic Bingham fluids and shear thinning and thickening H–B fluids in straight tubes of equilateral triangular cross-section. The asymptotic approach implemented does not involve expansions in terms of the Bingham number Bi. Thus, the solution method allows the mapping of the flow field of H–B fluids in tubes with axially symmetric but otherwise arbitrary contours for high Bi numbers. 2. Mathematical development The steady state field equations are, 𝒖 ⋅ ∇𝒖 = ∇ ⋅ 𝝈 , 𝝈 = −𝟏𝑝 + 𝝉 , ∇⋅𝒖=0 ∗ Corresponding author. E-mail addresses: [email protected] (M.F. Letelier), [email protected] (D.A. Siginer), [email protected] (E. González). https://doi.org/10.1016/j.ijnonlinmec.2023.104408 Received 5 January 2023; Received in revised form 16 March 2023; Accepted 3 April 2023 Available online 6 April 2023 0020-7462/© 2023 Elsevier Ltd. All rights reserved. M.F. Letelier, D.A. Siginer and E. González International Journal of Non-Linear Mechanics 153 (2023) 104408 We note that no assumptions were made up to this point on the shape of the cross-sectional contour shape of the tube that is no numerical value has been assigned to the shape parameter n up to this point. Thus, Eq. (9) governs at order 𝑂 (𝜀) of the asymptotic expansion the flow of H–B family of fluids in all tubes of axially symmetric cross-sectional contours represented by the mapping (2). The function 𝑓1 embedded in the asymptotic expansion (3) is computed through the solution of (9). Eq. (9) is also the governing equation for the flow of the power-law fluids for Bi=0, where (𝒖, 𝝈, 𝟏) and p represent the velocity vector, total stress tensor, the unit tensor and the pressure field, respectively. The H–B constitutive equation reads as, √ ( ) 𝜏𝑦 1 𝝉 = 𝜂II𝛼−1 𝑫 , II = 𝑫∶𝑫 (1) + 𝑫 𝑫 II𝑫 2 𝝉 and 𝜏𝑦 are the extra-stress tensor and the yield shear stress and II𝑫 and 𝜂 represent the second invariant of the rate of deformation tensor D and the viscosity, respectively. Herschel–Bulkley fluids can exhibit shear-thinning and shear thickening for values of the power index 𝛼 < 1 and 𝛼 > 1, respectively, in contrast to the Bingham viscoplastic fluids (𝛼 = 1) the viscosity 𝜂 of which is shear rate independent. The following scale factors are introduced together with the non-dimensional Bingham number 𝐵𝑖, 𝜂𝑤0 ∗ 𝜂𝑤0 ∗ 𝜏 , 𝑝= 𝑝 , 𝑡 = 𝑡0 𝑡∗ , 𝑟 = 𝑎𝑟∗ , 𝑤 = 𝑤0 𝑤∗ , 𝜏𝑟𝑧 = 𝑎 𝑟𝑧 𝑎 𝑎𝜏𝑦 𝜂𝑤0 𝜌𝑎2 , 𝑡0 = 𝐵𝑖 = , 𝜏0 = 𝜂𝑤0 𝜂 𝑎 𝑛 𝐺 (𝑟, 𝜃) = 1 − 𝑟 + 𝜀𝑟 sin 𝑛𝜃 (2) The velocity field is structured as, [ ( )] 𝑤 (𝑟, 𝜃) = 𝐺 (𝑟, 𝜃) 𝑓0 (𝑟) + 𝜀𝑓1 (𝑟, 𝜃) + 𝜀2 𝑓2 (𝑟, 𝜃) + 𝑂 𝜀3 (3) Substitution of (10) into (9) yields the following equation for 𝑅(𝑟) 𝛼𝑟 (2𝑟 − 𝐵𝑖) 𝑅′′ (𝑟) + [2 (2𝛼 − 1) 𝑟 − 𝛼𝐵𝑖] 𝑅′ (𝑟) − 18𝑅 (𝑟) = 0 𝛾 is a constant whose value is fixed by the regularity condition at 𝑟 = 1 as 𝛾= 𝑅 ⟶ 2 + [2 (2𝛼 − 1) 𝑟 − 𝛼𝐵𝑖] 𝜕 2 𝑤1 𝜕𝑤1 +2 =0 𝜕𝑟 𝜕𝜃 2 1−𝛼 𝛼 𝐵𝑖→0 𝑟−𝑎 Similarly, 𝑅(𝑟) collapses onto the 𝑂(𝜀) solution for the purely Bingham flow in the triangular cross-section when 𝛼 = 1. ) ( 9𝐵𝑖 2 𝐵𝑖3 9𝐵𝑖2 + 𝑟− 𝑟 + 𝑟3 (13) 𝑅 (𝑟) ⟶ − 𝛼→1 80 40 10 To shed more light on the method that underlies the approach to the solution implemented in this paper the solution for the Bingham fluids is developed in some detail as it is simpler than that of, but embedded in, the H–B fluids. In this case Eq. (11) is reduced to, 𝑟 (2𝑟 − 𝐵𝑖) 𝑅′′ (𝑟) + [2𝑟 − 𝐵𝑖] 𝑅′ (𝑟) − 18𝑅 (𝑟) = 0 (8) (14) To arrive at this equation starting from Eqs. (4) and (5) we set 𝛼 = 1 and substitute (6) in the Eqs. (4) and (5) to obtain 𝜏𝑟𝑧 = 𝐵𝑖− 1 − 𝑟2 Shear stress components can be linearized through a combination of (4)–(6) and substituted into (8). The resulting equation for 𝑤1 (𝑟, 𝜃) is 𝜕𝑟2 (2 − 𝐵𝑛)1∕𝛼 ) ( 2 2𝐹2,1 𝑎, 𝑏, 𝑐; 𝐵𝑖 We note that for 𝐵𝑖 = 0 the function 𝑅(𝑟) reduces to the 𝑂(𝜀) contribution to the purely power-law fluid solution, i.e. 𝑤0 (𝑟) 𝜕 2 𝑤1 (11) The solution to this equation is given by, ) ( 2𝑟 (12) 𝑅 (𝑟) = 𝛾𝐹2,1 𝑎, 𝑏, 𝑐; 𝐵𝑖 ( ) 2𝑟 where 𝐹2,1 𝑎, 𝑏, 𝑐; 𝐵𝑖 is the Hypergeometric function of the second kind with 1−𝛼−𝛥 1−𝛼+𝛥 ;𝑏 = − ;𝑐 = 1 𝑎=− 2𝛼 2𝛼 √ 𝛥 = (1 − 𝛼)2 + 36𝛼 Here 𝑤0 is the well-known velocity field of a H–B fluid in a circular tube, i.e. 𝛼 𝑤0 (𝑟) = 𝐾 [(𝑌 (1)) − 𝑌 (𝑟)] 𝐾= (7) 2 (1 + 𝛼) 1 + 𝛼 𝑌 (𝑟) = (2𝑟 − 𝐵𝑖)𝛽 𝛽= 𝛼 The equation of motion in terms of the extra-stress components and the modified pressure is given as, 𝛼𝑟 (2𝑟 − 𝐵𝑖) (10) 𝑤1 (𝑟, 𝜃) = 𝑅 (𝑟) sin 3𝜃 𝜏𝑟𝑧 and 𝜏𝜃𝑧 are the components of the extra-stress tensor, and 𝐵𝑖 and 𝛼 represent the dimensionless yield stress and the power index, respectively. In addition to the expansion (3) 𝑤(𝑟, 𝜃) is also expanded in an asymptotic series in terms of the small parameter 𝜀, ( ) 𝑤 = 𝑤0 (𝑟) + 𝜀𝑤1 (𝑟, 𝜃) + 𝜀2 𝑤2 (𝑟, 𝜃) + 𝑂 𝜀3 (6) 𝑓0 = + [2 (2𝛼 − 1) 𝑟] where 𝑅(𝑟) is a function of the radial coordinate 𝑟. If the mapping parameter n is assigned the value of 3 an equilateral triangle contour is obtained. The square is represented by 𝑛 = 4. We now proceed to seek as a specific example the solution of (9) for the description of the H–B fluid flow in a tube with equilateral triangular cross-section at order 𝑂 (𝜀) by prescribing the velocity field, where 𝑤(𝑟, 𝜃) is the axial velocity expressed in terms of the polar coordinates (𝑟, 𝜃), 𝜀 and 𝑛 represent geometrical mapping parameters and (𝑓0 , … , 𝑓𝑛 ) are unknown functions to be determined. The nonvanishing extra-stress components of a H–B fluid for parallel, 2-D flow, in dimensionless variables, are given in polar coordinates, ( ) ) ( 𝐵𝑖 𝜕𝑤 𝐵𝑖 1 𝜕𝑤 𝛼−1 𝜏𝑟𝑧 = II𝛼−1 + 𝜏 = II + (4) 𝜃𝑧 𝑫 𝑫 II𝑫 𝜕𝑟 II𝑫 𝑟 𝜕𝜃 1 [( ) ( ) ] 1 𝜕𝑤 2 2 𝜕𝑤 2 + II𝑫 = (5) 𝜕𝑟 𝑟 𝜕𝜃 𝜏𝑟𝑧 𝜕𝜏𝑟𝑧 1 𝜕𝜏𝜃𝑧 𝜕𝑝 + + =− =4 𝑟 𝜕𝑟 𝑟 𝜕𝜃 𝜕𝑧 The comparison of 𝑂(1) terms in (3) and (6) yields, 𝜕𝑟2 𝑤1 (𝑟, 𝜃) = 𝑅 (𝑟) sin 𝑛𝜃 where 𝑎, 𝑤0 , 𝑡0 and 𝜏0 are reference values for length, velocity, time and shear stress, respectively. The star in the subscript indicates dimensionless variables. For the sake of brevity, the star notation is omitted in the rest of the paper. However, all variables are dimensionless unless indicated otherwise. We consider parallel, laminar and isothermal viscoplastic flows in straight long tubes. The cross-sectional contour 𝐺 (𝑟, 𝜃) = 0 is determined through a continuous and one-to-one mapping 𝐺(𝑟, 𝜃) of the circular tube boundary of non-dimensional unit radius that satisfies the no-slip condition at the wall, 2 𝜕 2 𝑤1 𝜕𝑤1 𝜕 2 𝑤1 =0 +2 𝜕𝑟 𝜕𝜃 2 in the same tubes whose contours are in the same spectrum of axially symmetric tube contours represented by (2). The general solution of (9) for all contours represented by (2) is given by 2𝛼𝑟2 ( ) 𝜕𝑤0 𝜕𝑤 −𝜀 1 +𝑂 𝜀2 𝜕𝑟 𝜕𝑟 𝜏𝜃𝑧 = 𝜀 𝑟 ( ) 𝜕𝑤0 𝐵𝑖 − 𝜕𝑟 𝜕𝑤1 𝜕𝑟 𝜕𝑤0 𝜕𝑟 ( ) +𝑂 𝜀2 (15) The zeroth order velocity 𝑤0 in (6) is the well-known velocity of the Bingham fluid in round tubes except for the region of the plug zone, (9) 𝑤0 (𝑟) = 1 − 𝑟2 − 𝐵𝑖 (1 − 𝑟) 2 (16) M.F. Letelier, D.A. Siginer and E. González International Journal of Non-Linear Mechanics 153 (2023) 104408 which can also be obtained from the velocity field of the H–B fluid (7) as well by setting 𝛼 = 1. The continuously varying velocity expression (16) in the radial direction is not valid over the plug zone, the body of viscoplastic fluid around the axis that moves as a solid body. The location of the interface of the plug zone with the moving fluid at this order 𝑂(1) is determined by 𝑑 𝑤 (𝑟) = 0 𝑑𝑟 0 which sets the plug zone radius 𝑟0 and the plug velocity, problem. All regularization formulations approach the Bingham model as 𝜀 → 0 as in the Papanastasiou model given below, ( II ) ] [ ⎤ ⎡ − 𝑫 𝜏0 1 − 𝑒 𝜀 ⎥ ⎢ ⎥𝑫 𝝉 = 2 ⎢𝜂 + ⎥ ⎢ II𝑫 ⎥ ⎢ ⎦ ⎣ (17) Pending issues and difficulties with the simulation of viscoplastic flows are detailed in a relatively recent contribution, Denn and Bonn [8]. In the following, we present twelve figures representing the pressure gradient driven flow of viscoplastic Bingham and Herschel–Bulkley type of fluids together with power-law fluids in a straight tube of equilateral triangular cross-section. Each figure includes a plot of equal velocity lines and a corresponding velocity profile drawn along a symmetry axis. Plug zones in all the Figs. were determined using the necessary and sufficient condition (18), the equivalent of (17) at order 𝑂 (𝜀), ( ) 𝐵𝑖 𝐵𝑖2 , 𝑤0 𝑟0 = 1 − 𝐵𝑖 + = 𝐶𝑜𝑛𝑠𝑡. 2 4 Eq. (16) does not describe the flow field in the region 𝑟0 > 𝑟 ≥ 0. It shows a continuous variation in the plug zone ending up with an ( ) inverted spike and cusp [𝑤0 (0) = (1 − 𝐵𝑖) < 𝑤0 𝑟0 ] at the center ( ) 𝑟 = 0 smaller than 𝑤0 𝑟0 , which contradicts the physical evidence of the existence of the plug which moves as a solid body in the region 𝑟0 ≥ 𝑟 ≥ 0. However, (16) perfectly represents the velocity variation over the rest of the cross-section between the interface and the wall. A single velocity equation valid over the full circular cross-section is not available. However, the physical existence of the plug zone and the description (16) of the velocity of the freely flowing viscoplastic fluid between the interface and the wall are well established and accepted. Comparing 𝑂 (1) terms in the velocity fields defined in (3) and (6) and using (16) we determine the function 𝑓0 (𝑟) in the velocity expansion (3), 𝑟0 = 𝑓0 (𝑟) = 1 − 𝜕 𝑑 𝜕𝑤 = [𝑤0 (𝑟) + 𝜀𝑤1 (𝑟, 𝜃)] = −2𝑟 + 𝐵𝑖 + 𝜀𝑠𝑖𝑛3𝜃 𝑅 (𝑟) = 0 (18) 𝜕𝑟 𝜕𝑟 𝑑𝑟 The location 𝑟0 (𝜃) of the interface of the plug zone and the freely flowing viscoplastic fluid up to and including order 𝑂(𝜀) satisfies (18). 𝑅 (𝑟) in (18) is defined in (12). We begin with a comparison of the velocity field of a power-law shear-thinning fluid (𝛼 = 0.5, 𝐵𝑖 = 0.0) in Fig. 1 with the velocity field of a H–B fluid of the same shear-thinning capability as the powerlaw fluid in Fig. 1 but endowed with viscoplasticity in Fig. 2 (𝛼 = 0.5, 𝐵𝑖 = 0.5). We note the stagnation zones emerging in the corners of the cross-section in Fig. 2 and the much-reduced maximum velocity of the H–B fluid, about half of that of the power law fluid of the same shear-thinning capability in Fig. 1. The velocity distribution of the power-law fluid throughout the cross-section is smooth whereas the emergence of a cusp is noticeable at the maximum velocity point, the center of the triangle in the case of the shear-thinning viscoplastic fluid in Fig. 2. These anomalous features in the case of the shear-thinning viscoplastic fluid due to the interaction of the shear rate dependent viscosity (𝛼 < 1) and plasticity (𝐵𝑖 ≠ 0) become more pronounced as the value of the plasticity increases, Fig. 3& 4. The velocity profile along the symmetry axis shows that the velocity is reduced over much of the cross-section with a velocity spike close to the center and a cusp at the center. However, the later features are only mathematical artefacts of the velocity distribution ( ) 𝑤 = 𝑤0 (𝑟) + 𝜀𝑤1 (𝑟, 𝜃) + 𝑂 𝜀2 (19) 𝐵𝑖 1+𝑟 Substitution of (15) and (16) into (8) yields Eq. (14) the solution of which is given by (13). Comparison of 𝑂(𝜀) terms in (3) and (6) yields 𝑓1 (𝑟, 𝜃) in (3), 𝑓1 (𝑟, 𝜃) = 𝑅(𝑟) − 𝑟3 𝑓0 (𝑟) sin 3𝜃 1 − 𝑟2 The second order velocity 𝑤2 (𝑟, 𝜃) can be computed from ) ( 2 𝜕 2 𝑤2 𝜕 𝑤2 1 𝜕𝑤2 + = 2𝐾(𝑟, 𝜃) +2 𝑟(2𝑟 − 𝐵𝑖) 2 𝑟 𝜕𝑟 𝜕𝑟 𝜕𝜃 2 Where 𝐾(𝑟, 𝜃) is a complex function involving 𝑤0 (𝑟) and 𝑤1 (𝑟, 𝜃). Nevertheless, 𝑤2 (𝑟, 𝜃) can be computed analytically. But, the outcome does not change the values of the field variables qualitatively and are not included in the results reported in this paper. taken as valid and continuous over the entire cross-section whereas it is only valid over the section of the tube between the wall and the interface of the plug zone. That is (19) does not apply between the point 𝑟0 (𝜃) of the solution of (18) and the center 𝑟 = 0 between which points the velocity is uniform and equal to [ ] ( ) [ ] 𝑤 𝑟0 (𝜃) = 𝑤0 𝑟0 + 𝜀𝑤1 𝑟0 (𝜃) , 𝜃 = 𝐶𝑜𝑛𝑠𝑡. 3. Results and discussion Setting 𝐺(𝑟, 𝜃) = 0 and varying the values of the mapping parameters 𝜀 and 𝑛 yields a great variety of axially symmetric cross-sectional contours. In particular setting 𝜀 = 0.3849 and 𝑛 = 3 leads to an exact equilateral triangle. In Cartesian coordinates 𝐺 (𝑟, 𝜃) → 𝐺 (𝑥, 𝑦) reads as, ( √ )( √ ) 3 3 𝐺 (𝑥, 𝑦) = (3𝜀𝑦 − 1) 𝑥 + 1 + 𝑦 𝑥−1− 𝑦 3 3 We did not suppress the artefacts in Figs. 2–4 and chose to show them for clarity in one case. However, they are all suppressed in all the remaining Figs. and only the uniform physical velocity of the plug is shown. But, we stress again that this is a physical fact and not a mathematical outcome. We note that the mathematical artefacts in the plug region at order 𝑂(𝜀) are the opposite of those at order 𝑂(1). In the former the velocity spike and cusp are up corresponding to a sharp increase of velocity at the center and in the latter the opposite occurs the spike and cusp are in the opposite direction at the center corresponding to a sharp decrease in velocity. All these anomalies are suppressed in the remaining Figs. and only the physical uniform plug velocity is shown, except in Figs. 2–4 they are not suppressed for clarity. We note the increasing size of the stagnation zones in Figs. 3 and 4. The dynamics of the stagnation zones was explored by Letelier et al. [7]. They show analytically that the stagnation zone near the which is exactly the equation describing an equilateral triangle when 𝐺 (𝑥, 𝑦) = 0. We note that if 𝜀 = 0.35 instead of 𝜀 = 0.3849 one obtains an equilateral triangle with slightly rounded corners and if 𝜀 = 0.22 and 𝑛 = 4 we get a square with slightly rounded corners. It is also possible to represent exact square and rectangular contours with sharp corners using a slightly modified mapping. We note also that the regularization methods adopted in numerical studies to remove the discontinuity such as those of Bercovier and Engelman [4] and Papanastasiou [5] are not needed in analytical studies. The regularization methods, which involve a small parameter 𝜀, convert the flow field of the viscoplastic fluid into that of a purely viscous Newtonian fluid where the change in the viscosity is continuous over the whole field. Numerical computations rest on varying the small parameter 𝜀 to get convergence and recover the original discontinuous 3 M.F. Letelier, D.A. Siginer and E. González International Journal of Non-Linear Mechanics 153 (2023) 104408 Fig. 1. Velocity level lines and velocity profile along a symmetry axis for a shear-thinning power-law fluid: 𝛼 = 0.5, 𝐵𝑖 = 0.0. Fig. 2. Velocity level lines and velocity profile along a symmetry axis for a shear-thinning viscoplastic fluid: 𝛼 = 0.5, 𝐵𝑖 = 0.5. The small spike and cusp at the center are artefacts of assuming the velocity 𝑤(𝑟, 𝜃) is continuous over the full cross-section 𝑟𝑤 ≥ 𝑟 ≥ 0 and are only shown for clarity. 𝑤(𝑟, 𝜃) is not valid over 𝑟0 (𝜃) > 𝑟 ≥ 0 beyond the interface fluid. The plug in this Fig. moves at a constant dimensionless velocity 𝑤[𝑟0 (𝜃) , 𝜃] ≅ 0.55 where the locations of 𝑟0 along the 𝑟0 (𝜃) of the plug with the free flowing viscoplastic ( ) symmetry axis shown are at 𝑟0 ( 3𝜋 ) ≅ 0.5, 𝑟0 2 𝜋 2 ≅ 0.4. The fluid in the shaded areas is stagnant. fluid is smaller, which it should be, and the purely shear thinning fluid developed a velocity plateau around the center. Figs. 6–8 depict the velocity field of a purely viscoplastic fluid with a constant viscosity. These Figs. should be contrasted with Figs. 2–4, which do represent fluids with the same plasticity but with shear rate dependent viscosity. Two observations are important. Given a value of the plasticity coupled with shear rate dependent viscosity stagnation zones form at the corners but there are no plug zones centered around the axis moving as a solid body. However, for the same value of the plasticity coupled with constant viscosity there are no stagnation zones, but initially almost cylindrical plug zones form around the axis (𝛼 = 1.0, 𝐵𝑖 = 0.5) as in Fig. 6, growing increasingly larger and further deformed away from cylindrical shape with increasing plasticity as shown in Figs. 7 and 8 corners due to viscoplastic effects always have a convex boundary. The size of the stagnant zone depends on the Bingham number 𝐵𝑖 and the vertex angle that defines the separation between the walls, the bigger the Bingham number the larger the stagnation zone for the same vertex angle. Indeed, that is what we observe in Figs. 2–4. The shape of the boundary of the stagnant zone has been a source of controversy in the past as some investigators predicted a concave boundary Roquet and Saramito [9] and others a convex boundary Papanastasiou and Boudouvis [10] including the present authors Letelier et al. [4,5,11]. Letelier et al. [7] conclusively show that the boundary is indeed convex, and that is what we find in this investigation. Fig. 5 represents the velocity field of a Newtonian fluid (𝛼 = 1.0, 𝐵𝑖 = 0.0). Comparison with the velocity field of a purely shear thinning fluid in Fig. 1 shows that the maximum velocity of the Newtonian 4 M.F. Letelier, D.A. Siginer and E. González International Journal of Non-Linear Mechanics 153 (2023) 104408 Fig. 3. Velocity level lines and velocity profile along a symmetry axis for a shear-thinning viscoplastic fluid: 𝛼 = 0.5, 𝐵𝑖 = 1.0. The small spike and cusp at the center are artefacts of assuming the velocity 𝑤(𝑟, 𝜃) is continuous over the full cross-section 𝑟𝑤 ≥ 𝑟 ≥ 0 and are only shown for clarity. 𝑤(𝑟, 𝜃) is not valid over 𝑟0 (𝜃) > 𝑟 ≥ 0 beyond the interface fluid. The plug in this Fig. moves at a constant dimensionless velocity 𝑤[𝑟0 (𝜃) , 𝜃] ≅ 0.15 where the locations of 𝑟0 (𝜃) along the 𝑟0 (𝜃) of the plug with the free flowing viscoplastic ( ) symmetry axis shown are at 𝑟0 ( 3𝜋 ) ≅ 0.7, 𝑟0 2 𝜋 2 ≅ 0.6. The fluid in the shaded areas is stagnant. Fig. 4. Velocity level lines and velocity profile along a symmetry axis for a shear-thinning viscoplastic fluid: 𝛼 = 0.5, 𝐵𝑖 = 1.3. The small spike and cusp at the center are artefacts of assuming the velocity 𝑤(𝑟, 𝜃) is continuous over the full cross-section 𝑟𝑤 ≥ 𝑟 ≥ 0 and are only shown for clarity. 𝑤(𝑟, 𝜃) is not valid over 𝑟0 (𝜃) > 𝑟 ≥ 0 beyond the interface of the plug 𝑟0 (𝜃) with the free flowing viscoplastic ( ) fluid. The plug in this Fig. moves at a constant dimensionless velocity 𝑤[𝑟0 (𝜃) , 𝜃] ≅ 0.05 where the locations of 𝑟0 (𝜃) along the symmetry axis shown are at 𝑟0 ( 3𝜋 ) ≅ 0.8, 𝑟0 2 𝜋 2 ≅ 0.7. The fluid in the shaded areas is stagnant. with plasticities of 𝐵𝑖 = 1.0, 1.3, respectively. The velocity profile is smooth except for the discontinuity at the edge of the plug zone. Some comments on the apparent lack of stagnation zones in Figs. 6– 8 are called for. To begin with we note that the existence of the stagnation zones in the case of purely viscoplastic fluids was proven much earlier, Mosolov and Mjasnikov [12]. Next, we observe that the stagnation zone sizes in pressure gradient driven flow for the same opening angle at the corner and the same dimensionless pressure gradient in Letelier et al. [7] are very small even for the largest Bingham number 𝐵𝑖 = 0.2 used. Granted the problems studied in [7] and in this paper are different. Pressure gradient driven flow between two infinite flat plates that meet at an acute angle is studied in [7] whereas confined flow in a straight tube is investigated in the present paper. It can be argued that in confined pressure gradient driven flow the stagnation zone for the same opening angle, the same pressure gradient and Bingham number will indeed be much smaller than in the case of the flow between infinite plates meeting at an acute angle. Even for much larger 𝐵𝑖 =1.3 the stagnation zone is not detectable in Fig. 8. However, at this value of the 𝐵𝑖 the plug zone is already quite large and is covering half of the cross-section if not more. With further increases the 𝐵𝑖 may exceed the critical Bingham number 𝐵𝑖𝑐 for that particular cross-section and may end up choking the flow before the 5 M.F. Letelier, D.A. Siginer and E. González International Journal of Non-Linear Mechanics 153 (2023) 104408 Fig. 5. Velocity level lines and velocity profile along a symmetry axis for a Newtonian fluid: 𝛼 = 1.0, 𝐵𝑖 = 0.0. Fig. 6. Velocity level lines and velocity profile along a symmetry axis for a Bingham viscoplastic fluid: 𝛼 = 1.0, 𝐵𝑖 = 0.5. The shaded area moves at constant velocity as a rigid body (plug). appearance, however small, of a stagnation zone occurs. The critical Bingham number 𝐵𝑖𝑐 is not known for the viscoplastic fluid flow in a tube with equilateral triangular cross-section. However, Mosolov and Mjasnikov [13] proved using variational methods that flow stops in a square cross-sectional tube if 𝐵𝑖 exceeds 𝐵𝑖𝑐 , 𝐵𝑖 > 𝐵𝑖𝑐 and computed the value of the 𝐵𝑖𝑐 with no-slip on the walls of the square cross-section, 𝐵𝑖𝑐 = 4 √ 2+ = 1.0603 the boundary in a square cross-section was investigated numerically by ( ) Roquet and Saramito [14]. They determined that 𝐵𝑖𝑐 𝑠 depends on the degree of slip and introduced a slip parameter S. They find that ( ) 𝐵𝑖𝑐 𝑠 for S = 0.6 is twice that of the value (20) of 𝐵𝑖𝑐 with no-slip on the boundary in a square cross-sectional straight tube. However, our research in progress indicates that slip in non-circular cross-sectional tubes is not uniform on the boundary in contrast to circular tubes. The present authors are in the process of investigating analytically slip in non-circular cross-sections and results will be presented at a later date. We are also investigating at this time the critical Bingham number 𝐵𝑖𝑐 for the equilateral triangular cross-sectional pipes for different values of the power index 𝛼 in the Herschel–Bulkley constitutive structure. A different value of 𝐵𝑖𝑐 at which flow is choked off corresponds to each and every value of 𝛼. Fig. 9 represents the field of a purely shear thickening power law fluid with no plasticity (𝛼 = 1.5, 𝐵𝑖 = 0.0). Velocity peak is much (20) 𝜋 We expect the 𝐵𝑖𝑐 for the equilateral triangular cross-sectional tube to be larger than this value, which indeed seems to be the case as in Fig. 8 the flow has not been choked yet when 𝐵𝑖 = 1.3. It is also quite likely that the phase change viscoplastic material will slip at ( ) the boundary, which affects the value of the critical 𝐵𝑖𝑐 → 𝐵𝑖𝑐 𝑠 where s refers to slip. The flow of viscoplastic material with slip on 6 M.F. Letelier, D.A. Siginer and E. González International Journal of Non-Linear Mechanics 153 (2023) 104408 Fig. 7. Velocity level lines and velocity profile along a symmetry axis for a Bingham viscoplastic fluid: 𝛼 = 1.0, 𝐵𝑖 = 1.0. The shaded area moves at constant velocity as a rigid body (plug). Fig. 8. Velocity level lines and velocity profile along a symmetry axis for a Bingham viscoplastic fluid: 𝛼 = 1.0, 𝐵𝑖 = 1.3. The shaded area moves at constant velocity as a rigid body (plug). smaller, almost half, of that of the purely shear thinning power law fluid in Fig. 1 and does not display the semblance of a plateau centered around the axis as is the case in Fig. 1. The level lines of the velocity field for increasing values of the plasticity at the same value of the shear thickening power index 𝛼 are drawn in Figs. 10–12. Initially an almost cylindrical plug zone forms for the smallest value of the plasticity (𝐵𝑖 = 0.5) in Fig. 10. However, it grows with growing values of the plasticity assuming a pseudo triangular shape with rounded corners and bowed sides. Solid plug zone interface with the free flowing viscoplastic material is computed using the necessary and sufficient condition (18) in all the Figs. included in this paper. We note that for the same value of the Bi the plug zone in the case of the purely viscoplastic fluid is somewhat smaller than the plug zone in the case of the shear-thickening viscoplastic fluid, which leads to the tentative conclusion that flow could be chocked off earlier if the viscoplastic fluid is shear thickening rather than viscoplastic with constant viscosity. There appears to be no detectable stagnation zones in the range of the 𝐵𝑖 investigated as was the case for the purely viscoplastic fluid with constant viscosity. If the tentative conclusion is correct, that is if the 𝐵𝑖𝑐 for the purely viscoplastic fluid is larger than the 𝐵𝑖𝑐 for the viscoplastic shear thickening fluid, flow will be chocked off earlier in the case of the shear thickening fluid and by the same token the stagnation zones in the corners may not appear at all as the whole cross-section will be covered by the plug zone. The velocity field is continuous in the cross-section in all cases except for the discontinuity at the very edge of the plug. 4. Conclusion For any value of the yield stress 𝐵𝑖 and 𝛼 > 1 that is when the fluid is shear-thickening, velocity increases when compared to the case of 7 M.F. Letelier, D.A. Siginer and E. González International Journal of Non-Linear Mechanics 153 (2023) 104408 Fig. 9. Velocity level lines and velocity profile along a symmetry axis for a shear-thickening power-law fluid: 𝛼 = 1.5, 𝐵𝑖 = 0.0. Fig. 10. Velocity level lines and velocity profile along a symmetry axis for a shear-thickening viscoplastic fluid: 𝛼 = 1.5, 𝐵𝑖 = 0.5. The shaded area moves at constant velocity as a rigid body (plug). the Bingham fluid (𝛼 = 1) that is the plug moves at a larger velocity as compared to the case of purely viscoplastic Bingham material for the same value of the Bi. For shear-thinning H–B fluids, 𝛼 < 1, the opposite behavior occurs that is overall the velocity decreases when compared to Bingham fluid, that is for the same value of the Bi the plug moves at a slower pace. However, when the velocity is normalized with the average velocity, the opposite trend is true. These results are, in general, in agreement with the standard behavior of power-index fluids in circular pipes. A non-obvious result is that the plug-zone size and geometry are a function of the power index 𝛼. Its shape departs from a circle and becomes more complex as 𝛼 varies from unity. By the same token the Bingham number 𝐵𝑖𝑐 (𝛼) at which flow is chocked off is a function of the exponent 𝛼 in addition to its dependence on the cross-sectional contour. The stagnant zones also depend on 𝛼, they are only present in the flow field of the shear thinning viscoplastic fluids but are either absent or undetectably small in the flow field of the shear thickening viscoplastic or purely viscoplastic fluids in the range of the Bingham numbers Bi and the power index 𝛼 investigated. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability No data was used for the research described in the article. 8 M.F. Letelier, D.A. Siginer and E. González International Journal of Non-Linear Mechanics 153 (2023) 104408 Fig. 11. Velocity level lines and velocity profile along a symmetry axis for a shear-thickening viscoplastic fluid: 𝛼 = 1.5, 𝐵𝑖 = 1.0. The shaded area moves at constant velocity as a rigid body (plug). Fig. 12. Velocity level lines and velocity profile along a symmetry axis for a shear-thickening viscoplastic fluid: 𝛼 = 1.5, 𝐵𝑖 = 1.3. The shaded area moves at constant velocity as a rigid body (plug). References [9] N. Roquet, P. Saramito, An adaptive finite element method for viscoplastic flows in a square pipe with stick–slip at the wall, J. Non-Newton. Fluid Mech. 155 (3) (2008) 101–115. [10] T.C. Papanastasiou, A.G. Boudouvis, Flows of viscoplastic materials: Models and computations, Comput. Struct. 64 (1–4) (1997) 677–694. [11] M.F. Letelier, D.A. Siginer, J. Stöckle, C. Barrera, F. Godoy, C. Rosas, Bingham fluids: Deformation and energy dissipation in triangular cross section tube flow, Meccanica 53 (1–2) (2018) 161–173. [12] P.P. Mosolov, V.P. Mjasnikov, On stagnant flow regions of a visco-plastic medium in pipes, J. Appl. Math. Mech. 30 (1966) 705–719. [13] P.P. Mosolov, V.P. Mjasnikov, Variational methods in the theory of viscous-plastic medium, J. Appl. Math. Mech. 29 (1965) 468–492. [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-Newton. Fluid Mech. 155 (3) (2008) 101–115. [1] E.C. Bingham, An investigation of the laws of plastic flow, Bull. Bureau Stand. 13 (2016) 309–353. [2] E.C. Bingham, Fluidity and Plasticity, Mc Graw-Hill, New York, 1922. [3] W.H. Herschel, T. Bulkley, Measurement of consistency as applied to rubber-benzene solutions, Am. Soc. Test Proc. 26 (1926) 621–633. [4] M. Bercovier, M. Engelman, A finite-element method for incompressible non-Newtonian flows, J. Comput. Phys. 36 (1980) 313–326. [5] T.C. Papanastasiou, Flow of materials with yield stress, J. Rheol. 31 (1987) 385–404. [6] M.F. Letelier, D.A. Siginer, C.B. Hinojosa, On the physics of viscoplastic fluid flow in non-circular tubes, Int. J. Non-Linear Mech. 88 (2017) 1–10. [7] M.F. Letelier, D.A. Siginer, J. Stöckle, Stagnation zone near a corner in viscoplastic fluid flow, J. Fluids Eng. 144 (7) (2022) 071301, (9 pages). [8] M.M. Denn, D. Bonn, Issues in the flow of yield stress liquids, Rheol. Acta 50 (2011) 307–315. 9