Abstract
Identification of stagnant regions of viscoplastic fluid flows in production lines and equipment is of paramount importance owing to potential material degradation and process contamination. The present work introduces an assessment strategy to identify, classify and quantify unyielded regions with the objective of optimizing the flow conditions with the purpose of minimizing stagnant regions. Flow of Carbopol^{®} 980 in a Tbifurcation channel is adopted to illustrate the procedure. The rheological behavior of Carbopol^{®} 980 was simulated using the HerschelBulkley viscoplastic model regularized by Papanastasiou’s exponential approach. The analysis shows that three distinct types of stagnant unyielded regions take place in the bifurcation channel depending upon the Reynolds condition. Furthermore, the rheological characteristics of the fluid indicate the existence of an ideal Reynolds condition which allows the smallest flow stagnant area at the bifurcation zone.
Keywords:
Viscoplastic flow; HerschelBulkley fluid; Papanastasiou regularization; Tbifurcation channel
INTRODUCTION
Viscoplastic flows are found in distribution channels and equipment of a wide variety of production processes, such as in the cosmetic, pharmaceutical and food industries. Some examples of viscoplastic fluids commonly encountered are mayonnaise, ketchup, gelled products and toothpaste (Hammad, 2017Hammad, K.J. Inertial thermal convection in a suddenly expanding viscoplastic flow field. International Journal of Heat and Mass Transfer, 106, 829840 (2017). https://doi.org/10.1016/j.ijheatmasstransfer.2016.10.013
https://doi.org/10.1016/j.ijheatmasstran...
). The main characteristic of these materials is the existence of a yield stress, a critical value of stresses below which no flow takes place (Bird, et al., 1983Bird, R.B., Dai, G.C., Yarusso, B.J. The rheology and flow of viscoplasctic materials. Reviews in Chemical Engineering, 1, 170 (1983). https://doi.org/10.1515/revce19830102
https://doi.org/10.1515/revce19830102...
), so that both yielded and unyielded regions exist in the flow. The ideal viscoplastic models, such as Bingham plastic, HerschelBulkey and Casson, are discontinuous and numerical solutions for complex geometries require regularized models, such as the biviscosity equation (Tanner and Milthorpe, 1983Tanner, R.I., Milthorpe, J.F. Numerical simulation of the flow of fluids with yield stress. Proceedings of 3rd International Conference on Numerical Methods in Laminar and Turbulent Flow, Pineridge Press, Swansea, p.680690 (1983).) and models proposed by Bercovier and Engelman (1980Bercovier, M., Engelman, M. A finiteelement method for incompressible nonNewtonian flows. Journal of Computational Physics, 36, 313326 (1980). https://doi.org/10.1016/00219991(80)901631
https://doi.org/10.1016/00219991(80)901...
) and Papanastasiou (1987Papanastasiou, T.C. Flows of materials with yield. Journal of Rheology , 31, 385404 (1987). https://doi.org/10.1122/1.549926
https://doi.org/10.1122/1.549926...
). The latter was adopted in this study to investigate the hydrodynamic behavior of a HerschelBulkley fluid in a Tbifurcation. The Papanastasiou regularization uses an exponential modification now widely applied in numerical studies involving viscoplastic materials (Mitsoulis, 2007Mitsoulis, E. Flows of viscoplastic materials: models and computations. Rheology Reviews, 2007, 135178 (2007).).
In unyielded regions, the stress state falls below the yield stress threshold and the viscoplastic material behaves as a very highviscosity liquid, known also as plug region, which is transported by the flow. The size of the unyielded regions is associated with the Bingham number (Abdali et al., 1992Abdali, S.S., Mitsoulis, E., Markatos, N.C. Entry and exit flows of Bingham fluids. Journal of Rheology, 36, 389407 (1992). https://doi.org/10.1122/1.550350
https://doi.org/10.1122/1.550350...
; Mitsoulis et al., 1993Mitsoulis, E., Abdali, S.S., Markatos, N.C. Flow simulation of herschelbulkley fluids through extrusion dies. The Canadian Journal of Chemical Engineering, 71, 147160 (1993). https://doi.org/10.1002/cjce.5450710120
https://doi.org/10.1002/cjce.5450710120...
). Papanastasiou and Boudouvis (1997Papanastasiou, T.C., Boudouvis, A.G. Flows of viscoplastic materials: models and computations. Computers & Structures, 64, 677694 (1997). https://doi.org/10.1016/S00457949(96)001678
https://doi.org/10.1016/S00457949(96)00...
) evaluated (a) the nonshear regions in square duct flows by verifying the unyielded plug region (UPR) inside the channel (which is transported by a yielded flowing film (YFF)), and (b) an apparent unyielded region (AUR) located in the corners of the channel section. In cases of abrupt expansions and other complex geometries, regions of stagnation and deposition of material may occur in the channel (Mendes et. al., 2007Mendes, P.R.S., Naccache, M., Varges, P.R., Marchesini, F.H. Flow of viscoplastic liquids through axisymmetric expansionscontractions. Journal of NonNewtonian Fluid Mechanics , 142, 207217 (2007). https://doi.org/10.1016/j.jnnfm.2006.09.007
https://doi.org/10.1016/j.jnnfm.2006.09....
). Scott et al. (1998) also show that the existence of a yield stress reduces the size and strength of recirculation zones in channels with axisymmetric abrupt expansions and in a 180º planar bifurcation.
It is noteworthy that most recent studies featuring the simulation of viscoplastic flows have focused either on development of new numerical solution strategies or improvement of existing viscoplastic models to better fit experimental data. Such studies are important to establish the robustness of the discretization scheme, numerical method and material description; however, in most cases, actual rheological parameters and/or flow settings are rarely discussed. Contrastingly, this work aims to discuss important flow features of Carbopol^{®} 980 owing to its widespread use as a thickening and gelling agent in cosmetic and pharmaceutical industries.
Carbopol^{®} 980 is based on carboxyvinylic acids and presents transparent characteristics, being harmless and easy to prepare (Piau, 2007Piau, J.M. Carbopol gels: Elastoviscoplastic and slippery glasses made of individual swollen sponges: Mesoand macroscopic properties, constitutive equations and scaling laws. Journal of NonNewtonian Fluid Mechanics , 144, 129 (2007). https://doi.org/10.1016/j.jnnfm.2007.02.011
https://doi.org/10.1016/j.jnnfm.2007.02....
). Carbopol^{®} dispersions are found in many everyday products, from toothpastes to tile cleaners, and are also useful vehicles for functional ingredients (Roberts and Barnes, 2001Roberts, G.P., Barnes, H.A. New measurements of the flowcurves for Carbopol dispersions without slip artefacts. Rheologica Acta, 40, 499503 (2001). https://doi.org/10.1007/s003970100178
https://doi.org/10.1007/s003970100178...
). As a thickening solution used in the cosmetic industry, it has been largely applied in preparation of gels (Corrêa, 2005Corrêa, N.M., Camargo Jr., F.B., Ignácio, R.F., Leonardi, G.R. Avaliação do comportamento reológico de diferentes géis hidrofílicos. Revista Brasileira de Ciências Farmacêuticas, 41, 7378 (2005). https://doi.org/10.1590/S151693322005000100008
https://doi.org/10.1590/S15169332200500...
), with the purpose of emulsification, stabilization and rheological control (Kim, 2003Kim, J.Y., Song, J.Y., Lee, E.J., Park, S.K. Rheological properties and microstructures of Carbopol gel network system. Colloid and Polymer Science, 281, 614623 (2003). https://doi.org/10.1007/s0039600208087
https://doi.org/10.1007/s003960020808...
).
The onset of unyielded regions is typical of viscoplastic fluids owing to the existence of an initial yield stress. Use of such fluids (e.g. Carbopol^{®} 980) in cosmetic and pharmaceutical industries requires special care to avoid stagnant deposits inside the production lines and equipment. Stagnant regions can be formed at the solid walls and can cause material degradation and process contamination. The present work introduces an assessment strategy to identify, classify and quantify unyielded regions with the objective of optimizing the flow conditions in a Tbifurcation channel with the purpose of minimizing stagnant regions.
THEORETICAL FORMULATION
The physical model adopted in this work represents the steady state flow of a viscoplastic fluid through a plane channel with a Tbifurcation. Due to the yield stress of the rheological model, as the velocity field develops along the channel, a central region of small velocity variations (small shear rates) is formed, giving rise to a fluid behaviour known as plug flow. In regions close to the wall, the small shear rates lead to stagnant fluids. Therefore, the mathematical model, including the constitutive relation, must account for such flow features.
Governing equations: The general assumptions adopted to solve the problem are as follows: incompressible and laminar flow, steady state, negligible body forces, isothermal flow, and twodimensional geometry. The mass and momentum conservation equations in their conservative form are, respectively,
where τ_{ij} is the shear stress tensor, which, in the case of incompressible Newtonian fluids, is directly associated with the rate of strain tensor, D_{ij}, or else by the shear strain rate tensor, γ_{ij}, thus given by a constitutive equation,
where
is the rate of strain tensor and η is the apparent viscosity. Noticeably, the generalized Newtonian fluid equation was used to represent the dependence of viscosity on the equivalent strain rate (Bird et al., 1987Bird, R.B., Armstrong, R.C., Hassager, O. Dynamic of Polymeric Liquids. Fluid Mechanics, Vol. 1. 1987.), so that:
in which η(γ) is function of a scalar invariant of the strain rate tensor. For incompressible and shear flow cases (Bird et al., 1987Bird, R.B., Armstrong, R.C., Hassager, O. Dynamic of Polymeric Liquids. Fluid Mechanics, Vol. 1. 1987.):
is the equivalent shear rate.
Thus, the momentum equation for a viscoplastic fluid (based on the generalized Newton model) is
Constitutive modelling: Rheological studies have shown that Carbopol^{®} 980 has viscoplastic properties which can be described by the HerschelBulkley model (Rudert and Schwarze, 2009Rudert, A., Schwarze, R., Experimental and numerical investigation of a viscoplastic Carbopol gel injected into a prototype 3D mold cavity. Journal of NonNewtonian Fluid Mechanics , 161, 6068 (2009). https://doi.org/10.1016/j.jnnfm.2009.04.006
https://doi.org/10.1016/j.jnnfm.2009.04....
). The literature also shows that thixotropic characteristics of these gels can, in principle, be neglected (Coussot, 2014Coussot, P. Yield stress fluid flows: A review of experimental data. Journal of NonNewtonian Fluid Mechanics , 211, 3149 (2014). https://doi.org/10.1016/j.jnnfm.2014.05.006
https://doi.org/10.1016/j.jnnfm.2014.05....
). The general characteristic of ideal viscoplastic fluids is the existence of a yield stress threshold, τ_{0}, below which stresses are null. For instance, due to the yield stress, as the velocity field develops along a channel, a central region of small velocity variations is formed (small shear rates), giving rise to a fluid behavior known as plug flow. Although ideal models for viscoplastic fluids (e.g., Bingham, Casson and HerschelBulkley) do not present limitations to analytical solutions of simple problems (Bird et al., 1983Bird, R.B., Dai, G.C., Yarusso, B.J. The rheology and flow of viscoplasctic materials. Reviews in Chemical Engineering, 1, 170 (1983). https://doi.org/10.1515/revce19830102
https://doi.org/10.1515/revce19830102...
), some hindrances may arise for more complex geometries. Numerical approaches can overcome some problems; however, the nonlinearities of both governing equations and the constitutive relation impose significant difficulties even when numerical modeling is used. Burgos et al. (1999Burgos, G.R., Alexandrou, A.N., Entov, V. On the determination of yield surfaces in HerschelBulkley fluids. Journal of Rheology , 43, 463483 (1999). https://doi.org/10.1122/1.550992
https://doi.org/10.1122/1.550992...
) highlight the fact that the discontinuity of the constitutive relation constitutes an inherent obstacle to numerical approximations.
The difficulties arise from the necessity to map the transitional surface between the shear flow  noshear flow condition established by the magnitude of the local yield stress. In its classical form, the apparent viscosity becomes unbounded due to the presence of γ in the denominator of the viscoplastic equations. Furthermore, even though the velocity field is calculated, the shape and location of the transition region are unknown. The discontinuity of the constitutive relation associated with the yield stress leads to high values of the viscosity function for small shear rates, which causes in some cases η(γ)→∞, with consequent numerical instabilities (Min et al., 1997Min, T., Choi, H.G., Yoo, J.Y., Choi, H. Laminar convective heat transfer of a Bingham plastic in a circular pipeII: Numerical approach hydrodynamically developing flow and simultaneously developing flow. International Journal of Heat and Mass Transfer , 40, 36893701 (1997). https://doi.org/10.1016/S00179310(97)000045
https://doi.org/10.1016/S00179310(97)00...
).
Papanastasiou (1987Papanastasiou, T.C. Flows of materials with yield. Journal of Rheology , 31, 385404 (1987). https://doi.org/10.1122/1.549926
https://doi.org/10.1122/1.549926...
) proposed an exponential regularization for the ideal Bingham model by introducing a parameter “m” that controls an exponential increase of the stresses. If its value is sufficiently high, the regularized model approximates Bingham’s ideal fluid behavior. The regularized Bingham model can be expressed as:
where η_{0} is a constant viscosity parameter and τ_{0} is the yield stress. The Papanastasiou regularization applied to the HerschelBulkley viscoplastic model provides:
in which K and n are the consistency and power law indices, respectively.
Numerical methodology: The main focus of the present work is to study important aspects of the viscoplastic flow inside channels, i.e., possible formation of stagnant regions that may contaminate the production line in the cosmetic, pharmaceutical or food processing industries. Therefore, the numerical scheme chosen was classical and only a brief discussion is presented. The numerical simulations were performed using the ANSYS FLUENT^{®} commercial software. The program uses the finite volume method (FVM) to discretize the conservation equations based on a cellcentered formulation and a second order scheme for discretization of the viscous terms. The SIMPLE method was adopted to approach the pressurevelocity coupling. The oddeven decoupling was resolved using the RhieChow interpolation scheme. The viscosity function was determined using the regularized HerschelBulkeyPapanastasiou model implemented via UDF (User Defined Function). The solutions of the linear systems were considered converged when the normalized residue reached R^{φ} < 10^{7}.
RESULTS AND DISCUSSION
Verification and validation of the numerical solution: This section presents the validation of the numerical method. The numerical solution of a viscoplastic fluid flow in a twodimensional plane channel using ANSYS FLUENT^{®} (based on the finite volume method) was compared against solutions (i) obtained using an inhouse code developed by the authors based on finite differences (FDM) and (ii) reported by Boualit et al. (2011Boualit, A., Zeraibi, N., Boualit, S., Amoura, M. Thermal development of the laminar flow of a Bingham fluid between two plane plates with viscous dissipation. International Journal of Thermal Sciences, 50, 3643 (2011). https://doi.org/10.1016/j.ijthermalsci.2010.09.005
https://doi.org/10.1016/j.ijthermalsci.2...
) determined using the finite element method. The flow conditions, geometry and rheological parameters follow the analysis performed by Boualit et al. (2011), who adopted the BinghamPapanastasiou model to compute the apparent viscosity of the viscoplastic fluid.
Figure 1 shows the geometry of the problem. The parameters and boundary conditions for the hydrodynamic analysis are: uniform axial velocity in the channel inlet (Γ_{1}: u^{*} = 1, v^{*} = 0); nonslip condition at the channel walls (Γ_{2}: u^{*} = 0, v^{*} = 0); fullydeveloped velocity at the channel outlet (Γ_{3}: du^{*}/dx^{*} = 0, v^{*} = 0); symmetry condition at the center of the channel (Γ_{4}: du^{*}/dy^{*} = 0, v^{*} = 0); null pressure at the channel outlet (Γ_{3}), and null normal pressure gradient at the channel walls (Γ_{2}).
Validation geometry  plane channel (Boualit et al., 2011Boualit, A., Zeraibi, N., Boualit, S., Amoura, M. Thermal development of the laminar flow of a Bingham fluid between two plane plates with viscous dissipation. International Journal of Thermal Sciences, 50, 3643 (2011). https://doi.org/10.1016/j.ijthermalsci.2010.09.005
https://doi.org/10.1016/j.ijthermalsci.2... ).
In order to guarantee the fullydeveloped velocity distribution towards the exit, a channel length L = 80H was considered. A mesh containing 300 x 30 elements was also used in the simulations. A comparative assessment of the axial velocities for fullydeveloped flow for different Bingham numbers, Bn, against the results reported by Boualit et al (2011Boualit, A., Zeraibi, N., Boualit, S., Amoura, M. Thermal development of the laminar flow of a Bingham fluid between two plane plates with viscous dissipation. International Journal of Thermal Sciences, 50, 3643 (2011). https://doi.org/10.1016/j.ijthermalsci.2010.09.005
https://doi.org/10.1016/j.ijthermalsci.2...
) and those obtained using the FDM, as shown in Figure 2, indicates that the present numerical methodologies and implementation of the regularized viscosity function are able to recover the reference results with acceptable accuracy. Furthermore, it can be observed that higher values of the yield stress, τ_{0} (which is a rheological characteristic of the fluid and associated with Bn), cause an increase of a portion of the fluid moving with uniform velocity (UPR).
Table 1 shows the friction factor, fRe, for Bingham numbers ranging from Bn = 0 to Bn = 6.5. The results are compared with the numerical simulations from Boualit et al. (2011Boualit, A., Zeraibi, N., Boualit, S., Amoura, M. Thermal development of the laminar flow of a Bingham fluid between two plane plates with viscous dissipation. International Journal of Thermal Sciences, 50, 3643 (2011). https://doi.org/10.1016/j.ijthermalsci.2010.09.005
https://doi.org/10.1016/j.ijthermalsci.2...
) and the analytical study carried out by Quaresma and Macêdo (1998Quaresma, J.N.N., Macêdo, E.N. Integral transform solution for the forced convection of HerschelBulkley fluids in circular tubes and parallelplates ducts. Brazilian Journal of Chemical Engineering, 15, 7779 (1998). https://doi.org/10.1590/S010466321998000100008
https://doi.org/10.1590/S01046632199800...
). The maximum difference between the present finite volume approximation and the analytical results is 2.37%. The differences can be explained by the use of the Papanastasiou regularization model in comparison to the ideal Bingham model adopted by Quaresma and Macêdo (1998).
The length necessary to reach the fully developed condition was evaluated considering higher Reynolds number conditions. Figure 3 illustrates the development of the axial velocity in the center of the channel for a Bingham number Bn = 1. The inlet length is proportional to the Reynolds number and the results agree with the values presented by Boualit et al. (2011Boualit, A., Zeraibi, N., Boualit, S., Amoura, M. Thermal development of the laminar flow of a Bingham fluid between two plane plates with viscous dissipation. International Journal of Thermal Sciences, 50, 3643 (2011). https://doi.org/10.1016/j.ijthermalsci.2010.09.005
https://doi.org/10.1016/j.ijthermalsci.2...
).
Flow in a Tbifurcation channel  mesh refinement study: The distribution process of a viscoplastic material is simulated considering a Tbifurcation channel, as depicted in Figure 4. The Reynolds number (Re ≤ 50) and the channel length (L = 20H) used in the simulations prompt fully developed flow upstream from the bifurcation zone at the channel outlets. The half thickness of the inlet and outlet channels is H = 0.015 m.
The boundary conditions (see Figure 4) are: uniform inlet velocity (Γ_{1}: u = u_{c}, v = 0); fullydeveloped velocity distribution at the channel outlets (Γ_{2} and Γ_{3}: u = 0, dv/dy = 0); nonslip condition at the channel walls (Γ_{4}: u = 0, v = 0); null pressure at the channel exit (Γ_{2} and Γ_{3}: p = 0). The simulations were performed considering the rheological properties of Carbopol^{®} 980 according to experimental data and curve fitting reported by Rudert and Schwarze (2009Rudert, A., Schwarze, R., Experimental and numerical investigation of a viscoplastic Carbopol gel injected into a prototype 3D mold cavity. Journal of NonNewtonian Fluid Mechanics , 161, 6068 (2009). https://doi.org/10.1016/j.jnnfm.2009.04.006
https://doi.org/10.1016/j.jnnfm.2009.04....
). Table 2 presents the parameters adopted in the HerschelBulkleyPapanastasiou model for the Carbopol^{®} 980 viscoplastic fluid.
Parameters of the HerschelBulkleyPapanastasiou model for Carbopol^{®} 980 (Rudert and Schwarze, 2009Rudert, A., Schwarze, R., Experimental and numerical investigation of a viscoplastic Carbopol gel injected into a prototype 3D mold cavity. Journal of NonNewtonian Fluid Mechanics , 161, 6068 (2009). https://doi.org/10.1016/j.jnnfm.2009.04.006
https://doi.org/10.1016/j.jnnfm.2009.04.... ).
The value of the regularization parameter m was defined according to the works of Burgos and Alexandrou (1999Burgos, G.R., Alexandrou, A.N., Entov, V. On the determination of yield surfaces in HerschelBulkley fluids. Journal of Rheology , 43, 463483 (1999). https://doi.org/10.1122/1.550992
https://doi.org/10.1122/1.550992...
), Burgos et al. (1999), Alexandrou et al. (2001) and Boualit et al. (2011Boualit, A., Zeraibi, N., Boualit, S., Amoura, M. Thermal development of the laminar flow of a Bingham fluid between two plane plates with viscous dissipation. International Journal of Thermal Sciences, 50, 3643 (2011). https://doi.org/10.1016/j.ijthermalsci.2010.09.005
https://doi.org/10.1016/j.ijthermalsci.2...
). The Reynolds number for the HerschelBulkley viscoplastic fluid can be derived as:
whereas the corresponding relationship between the Bingham and Reynolds numbers is
Assessment of mesh dependency was evaluated by the dimensionless pressure drop, Δp^{*}, and fullydeveloped dimensionless velocities u^{*} (x = 15H) and v^{*} (y = ?20Η) at the center of the channel for a Bingham number Bn = 1. A structured mesh with refinement at the T junction, inlet and outlet regions was used, as illustrated in Figure 5. The mesh size used in this investigation comprised 850, 5400, 17100, 28125, 37800 and 53350 elements.
The values of the pressure drop and maximum flow velocities in the inlet and exit predefined locations in relation to the number of elements are shown in Figure 6. The simulations show that the mesh with 28,125 elements presents a relative difference of the pressure drop smaller than 0.01% when compared to the more refined case. The corresponding relative differences for velocities u^{*} (x = 15H) and v^{*} (y = ?20Η) in the center line of the channel are 0.07% and 0.12%, respectively. Figure 7 shows the evolution of the fullydeveloped axial velocity across the inlet channel at x = 15H. Assessment of the effects of mesh refinement was also performed for other inlet velocities with similar results. Therefore, based upon the results of the relative errors for the pressure drop and velocities in the fullydeveloped regions, the mesh containing 28,125 elements was selected for the study.
Assessment of mesh dependency: (a) maximum velocities at the channel center line at predefined locations and (b) pressure drop.
Analysis of the flow in a channel with Tbifurcation: Identification of unyielded flows, especially stagnant regions, is of prime importance in both (a) production line / equipment design and (b) establishment of corresponding flow conditions. This section introduces an objective strategy to identify, classify and quantify apparent unyielded regions (AURs) aiming at optimization of the flow conditions for the target geometry. The study is performed for Carbopol^{®} 980 (Table 2) viscoplastic fluid flow in the Tbifurcation channel depicted in Figure 5, for Reynolds numbers ranging from 0.1 to 50 (i.e. inlet velocities u_{c} = 0.02014 to 1.269 m/s).
In order to identify unyielded regions, the local shear stress is compared to the yield stress. Initially, the shear strain rate field is determined from velocities, followed by computation of the local shear stress (Equation 12) based on a similar strategy adopted in the work of Zisis and Mitsoulis (2002Zisis, T., Mitsoulis, E. Viscoplastic flow around a cylinder kept between parallel plates. Journal of NonNewtonian Fluid Mechanics , 105, 120 (2002). https://doi.org/10.1016/S03770257(02)000253
https://doi.org/10.1016/S03770257(02)00...
). Therefore, by using Equations (5) and (6), the magnitude of the local shear stress is
so that, if τ < τ_{0}, no shear flow takes place. Computation of the local yield stress was also implemented in the commercial code via UDF.
Identification of the type of unyielded regions is relevant when setting the flow conditions for a particular geometry. This work classifies the unyielded regions into four distinct types: (i) “AURs type A” are stagnant regions which arise due to flow split at solid walls; (ii) “AURs type B” are stagnant regions attached to solid walls which are formed in areas of low pressure caused by changes in flow direction; (iii) “AUR type C” are unyielded regions formed in recirculation zones close to solid walls; (iv) It is relevant to mention that plug flow takes place close to the channel centerline; notwithstanding, in this case, the unyielded material is carried by the flow and dissipated when the velocity gradient (strain rate) increases.
Figures 8, 9 and 10 show the development of the streamlines and AURs with respect to the Reynolds number. For visualization purposes, the range of the dimensionless shear stress field, τ^{*} = τ/τ_{0}, is clipped at τ^{*} < 1, so that the white color represents regions with no apparent shear stress. For a Reynolds number Re = 1 (Figure 8), the low velocities favor onset of the “AUR type A” at the stagnation point located at the root of the Tbifurcation. The simulations show that the size of the “AUR type A” decreases as the Reynolds number increases. Figure 9 indicates two “AURs type B” attached to the walls next to the corners for a Reynolds number Re = 30. For larger Reynolds numbers the “AURs type C” are formed at the recirculation zones close to the corner, as illustrated in Figure 10 for Re = 50.
Streamlines and shear stress field for Re = 0.1: the AUR type A is located at the stagnation point.
Streamlines and shear stress field for Re = 50: AUR type C takes place in the recirculation zone.
It is clear by superimposing the streamlines on the dimensionless shear stress field, especially in the case of Re = 50, that there are recirculation regions where τ^{*} < 1. The apparent contradiction (prediction of simultaneous recirculation flows and unyielded regions) is due to the use of the regularized constitutive model for viscosity. The continuity of the viscoplastic HerschelBulkleyPapanastasiou equation (Equation 9) for any magnitude of shear rates allows the determination of velocity fields in regions where the shear stress magnitude is smaller than the yield stress. The regularization is also related to the decentralization between the vortices and the nonshear region apparently in recirculation. These results agree with other works available in the literature, such as Alexandrou et al. (2001Alexandrou, A.N., Mcgilvreay, T.M., Burgos, G. Steady HerschelBulkley fluid flow in threedimensional expansions. Journal of NonNewtonian Fluid Mechanics, 100, 7796 (2001). https://doi.org/10.1016/S03770257(01)001276
https://doi.org/10.1016/S03770257(01)00...
), who also found similar results by evaluating a HerschelBulkley fluid in abrupt expansions and pondered that, qualitatively, the observed physics of the recirculation regions is more consistent with those of generalized Newtonian fluids than with ideal HerschelBulkley fluids. The parameters of regularized models simply adjust the magnitude of local stresses and strain rates to estimate ideal conditions, but can never predict absolute nonshear regions. Notwithstanding, the results from regularized functions are well accepted for analysis of the hydrodynamic behavior of viscoplastic fluids, making it possible to determine process conditions which provide a minimum amount of stagnant fluid in the distribution line.
Quantification of the AUR sizes for all different types is based on the dimensionless relative area, Ω^{(i)} _{AUR}, where τ^{*} < 1, as:
where i = A, B, C indicates the corresponding AUR type and respective areas A^{(i)} _{AUR}.
Figure 11 shows the composition of the apparent unyielded regions type A, B and C, quantified by the dimensionless area, Ω_{AUR}, for flow conditions represented by the Reynolds number. The simulations indicate that small Reynolds numbers (low velocities or narrow channels for Carbopol^{®} 980) are associated with larger regions of stagnant fluid at the root of the bifurcation channel (AUR type A). As the Reynolds number increases, velocity gradients also increase in such regions, causing the AUR type A to decrease and, eventually, vanish. On the other hand, higher Reynolds numbers cause the regions of low pressure near the corners to increase, leading to formation of AUR type B. This is the only type of stagnant region for Reynolds numbers within the interval 20 ≤ Re ≤ 40. It is interesting to note that, for Reynolds numbers Re > 40, the flow dynamics give rise to a detached unyielded region near each corner (AUR type C), causing the stagnant regions at the solid walls (AUR type B) to decrease (see also Figures 9 and 10).
The ideal flow condition was determined for a Reynolds number Re = 15, which provides the smallest stagnant area for the present Tbifurcation channel. This Re constitutes also the transition between formation of stagnant regions at the root of the bifurcation channel and at the corner of the exit channels. Therefore, for the present geometrical configuration, the flow conditions for Carbopol^{®} 980 corresponding to Re = 15 (inlet velocity u_{c} = 0.568 m/s) are highly recommended in order to minimize the chance of material degradation and process contamination.
CONCLUDING REMARKS
The flow characteristics of a viscoplastic fluid in a Ttype bifurcation were evaluated by numerical solution. The finite volume method was applied using the ANSYS FLUENT^{®} package associated with Papanastasiou’s regularization for HerschelBulkley fluids implemented via UDF. The numerical approach was verified by comparing hydrodynamic analytical and other numerical results for fluid flow between parallel plates. The rheological model implemented in the aforementioned commercial software was able to recover the expected hydrodynamic behavior with acceptable accuracy.
The numerical study of a viscoplastic flow was performed for the Carbopol^{®} 980 fluid, which is a thickener agent widely used in cosmetic and pharmaceutical industries. The flow condition was assessed for Reynolds numbers 0.1 ? Re ? 50 based on streamlines, shear stress fields and apparent unyielded regions (AURs).
This work also proposed an assessment strategy to identify, classify and quantify unyielded regions with the objective of optimizing the flow conditions with the purpose of minimizing stagnant regions. Based upon the procedure, the following conclusions can be ascertained:
In the inlet channel, upstream of the bifurcation zone, there is formation of an unyielded region (plug flow) where the fully developed condition is attained. In this region, the local shear stresses are smaller than the yield stress of the viscoplastic material. The unyielded material carried by the flow is dissipated when the velocity gradient increases at the bifurcation junction.
Stagnant regions are found at the root of the bifurcation where the flow splits between the exit channels. The total area of this type of unyielded material decreases when the Reynolds number (i.e. inlet velocities) increases.
Stagnant fluid attached to the exit walls and at recirculation regions near the bifurcation corner is found only for higher Reynolds numbers.
A Reynolds number Re = 15 prompted the recommended condition for Carbopol^{®} 980 flow in the Tbifurcation channel addressed in this work by providing the smallest fraction of all types of unyielded stagnant regions.
ACKNOWLEDGMENTS
The first author acknowledges CAPES (Brazilian Agency for Improvement of Higher Level Education Personnel) for the scholarship, Finance Code 001.
REFERENCES
 Abdali, S.S., Mitsoulis, E., Markatos, N.C. Entry and exit flows of Bingham fluids. Journal of Rheology, 36, 389407 (1992). https://doi.org/10.1122/1.550350
» https://doi.org/10.1122/1.550350  Alexandrou, A.N., Mcgilvreay, T.M., Burgos, G. Steady HerschelBulkley fluid flow in threedimensional expansions. Journal of NonNewtonian Fluid Mechanics, 100, 7796 (2001). https://doi.org/10.1016/S03770257(01)001276
» https://doi.org/10.1016/S03770257(01)001276  Bercovier, M., Engelman, M. A finiteelement method for incompressible nonNewtonian flows. Journal of Computational Physics, 36, 313326 (1980). https://doi.org/10.1016/00219991(80)901631
» https://doi.org/10.1016/00219991(80)901631  Bird, R.B., Armstrong, R.C., Hassager, O. Dynamic of Polymeric Liquids. Fluid Mechanics, Vol. 1. 1987.
 Bird, R.B., Dai, G.C., Yarusso, B.J. The rheology and flow of viscoplasctic materials. Reviews in Chemical Engineering, 1, 170 (1983). https://doi.org/10.1515/revce19830102
» https://doi.org/10.1515/revce19830102  Boualit, A., Zeraibi, N., Boualit, S., Amoura, M. Thermal development of the laminar flow of a Bingham fluid between two plane plates with viscous dissipation. International Journal of Thermal Sciences, 50, 3643 (2011). https://doi.org/10.1016/j.ijthermalsci.2010.09.005
» https://doi.org/10.1016/j.ijthermalsci.2010.09.005  Burgos, G.R., Alexandrou, A.N. Flow development of HerschelBulkley fluids in a sudden threedimensional square expansion. Journal of Rheology , 43, 485498 (1999). https://doi.org/10.1122/1.550993
» https://doi.org/10.1122/1.550993  Burgos, G.R., Alexandrou, A.N., Entov, V. On the determination of yield surfaces in HerschelBulkley fluids. Journal of Rheology , 43, 463483 (1999). https://doi.org/10.1122/1.550992
» https://doi.org/10.1122/1.550992  Corrêa, N.M., Camargo Jr., F.B., Ignácio, R.F., Leonardi, G.R. Avaliação do comportamento reológico de diferentes géis hidrofílicos. Revista Brasileira de Ciências Farmacêuticas, 41, 7378 (2005). https://doi.org/10.1590/S151693322005000100008
» https://doi.org/10.1590/S151693322005000100008  Coussot, P. Yield stress fluid flows: A review of experimental data. Journal of NonNewtonian Fluid Mechanics , 211, 3149 (2014). https://doi.org/10.1016/j.jnnfm.2014.05.006
» https://doi.org/10.1016/j.jnnfm.2014.05.006  Hammad, K.J. Inertial thermal convection in a suddenly expanding viscoplastic flow field. International Journal of Heat and Mass Transfer, 106, 829840 (2017). https://doi.org/10.1016/j.ijheatmasstransfer.2016.10.013
» https://doi.org/10.1016/j.ijheatmasstransfer.2016.10.013  Kim, J.Y., Song, J.Y., Lee, E.J., Park, S.K. Rheological properties and microstructures of Carbopol gel network system. Colloid and Polymer Science, 281, 614623 (2003). https://doi.org/10.1007/s0039600208087
» https://doi.org/10.1007/s0039600208087  Mendes, P.R.S., Naccache, M., Varges, P.R., Marchesini, F.H. Flow of viscoplastic liquids through axisymmetric expansionscontractions. Journal of NonNewtonian Fluid Mechanics , 142, 207217 (2007). https://doi.org/10.1016/j.jnnfm.2006.09.007
» https://doi.org/10.1016/j.jnnfm.2006.09.007  Min, T., Choi, H.G., Yoo, J.Y., Choi, H. Laminar convective heat transfer of a Bingham plastic in a circular pipeII: Numerical approach hydrodynamically developing flow and simultaneously developing flow. International Journal of Heat and Mass Transfer , 40, 36893701 (1997). https://doi.org/10.1016/S00179310(97)000045
» https://doi.org/10.1016/S00179310(97)000045  Mitsoulis, E. Flows of viscoplastic materials: models and computations. Rheology Reviews, 2007, 135178 (2007).
 Mitsoulis, E., Abdali, S.S., Markatos, N.C. Flow simulation of herschelbulkley fluids through extrusion dies. The Canadian Journal of Chemical Engineering, 71, 147160 (1993). https://doi.org/10.1002/cjce.5450710120
» https://doi.org/10.1002/cjce.5450710120  Papanastasiou, T.C. Flows of materials with yield. Journal of Rheology , 31, 385404 (1987). https://doi.org/10.1122/1.549926
» https://doi.org/10.1122/1.549926  Papanastasiou, T.C., Boudouvis, A.G. Flows of viscoplastic materials: models and computations. Computers & Structures, 64, 677694 (1997). https://doi.org/10.1016/S00457949(96)001678
» https://doi.org/10.1016/S00457949(96)001678  Piau, J.M. Carbopol gels: Elastoviscoplastic and slippery glasses made of individual swollen sponges: Mesoand macroscopic properties, constitutive equations and scaling laws. Journal of NonNewtonian Fluid Mechanics , 144, 129 (2007). https://doi.org/10.1016/j.jnnfm.2007.02.011
» https://doi.org/10.1016/j.jnnfm.2007.02.011  Quaresma, J.N.N., Macêdo, E.N. Integral transform solution for the forced convection of HerschelBulkley fluids in circular tubes and parallelplates ducts. Brazilian Journal of Chemical Engineering, 15, 7779 (1998). https://doi.org/10.1590/S010466321998000100008
» https://doi.org/10.1590/S010466321998000100008  Roberts, G.P., Barnes, H.A. New measurements of the flowcurves for Carbopol dispersions without slip artefacts. Rheologica Acta, 40, 499503 (2001). https://doi.org/10.1007/s003970100178
» https://doi.org/10.1007/s003970100178  Rudert, A., Schwarze, R., Experimental and numerical investigation of a viscoplastic Carbopol gel injected into a prototype 3D mold cavity. Journal of NonNewtonian Fluid Mechanics , 161, 6068 (2009). https://doi.org/10.1016/j.jnnfm.2009.04.006
» https://doi.org/10.1016/j.jnnfm.2009.04.006  Scott, P.S., Mirza, F., Vlachopoulos, J. Finiteelement simulation of laminar viscoplastic flows with regions of recirculation. Journal of Rheology , 32, 387400 (1988). https://doi.org/10.1122/1.549976
» https://doi.org/10.1122/1.549976  Tanner, R.I., Milthorpe, J.F. Numerical simulation of the flow of fluids with yield stress. Proceedings of 3^{rd} International Conference on Numerical Methods in Laminar and Turbulent Flow, Pineridge Press, Swansea, p.680690 (1983).
 Zisis, T., Mitsoulis, E. Viscoplastic flow around a cylinder kept between parallel plates. Journal of NonNewtonian Fluid Mechanics , 105, 120 (2002). https://doi.org/10.1016/S03770257(02)000253
» https://doi.org/10.1016/S03770257(02)000253
NOMENCLATURE
 AUR Apparent unyielded region
 A_{AUR} Apparent unyielded area [m^{2}]
 Bn Bingham number: Bn = (τ_{0}2H)/(ηu_{c})
 D_{ij} Rate of strain tensor [1/s]
 fRe Dimensionless friction factor
 H Characteristic inlet channel width [m]
 H_{E1} Exit channel width 1 [m]
 H_{E2} Exit channel width 2 [m]
 K Consistency index [Pa . s^{n}]
 L Channel length [m]
 m Papanastasiou’s regularization parameter [1/s]
 n Behavior (or powerlaw) index
 p Static pressure [Pa]
 p^{*} Dimensionless static pressure: p^{*} = (pH^{n})/(Ku_{c} ^{n})
 Δp^{*} Dimensionless static average pressure
 Re Reynolds number: Re = (ρu_{c} ^{2n}H_{n})/K
 R^{ϕ} Normalized residue
 t Time [s]
 u Horizontal component of the velocity [m/s]
 u^{*} Dimensionless horizontal component of the velocity: u^{*} = u/u_{c}
 u_{ij} Characteristic velocity [m/s]
 u_{i}; u_{j} Cartesian components of the velocity [m/s]
 UDF User Defined Function
 UPR Unyielded plug region
 v Vertical component of the velocity [m/s]
 v^{*} Dimensionless vertical component of the velocity: v^{*} = v/u_{c}
 x Horizontal coordinate [m]
 x^{*} Dimensionless horizontal coordinate: x^{*} = x/H
 x_{i}; x_{j} Cartesian components [m]
 y Vertical coordinate [m]
 y^{*} Dimensionless vertical coordinate: y^{*} = y/H
 YFF Yielded flowing film

Greek letters
 η Aapparent viscosity [Pa . s]
 γ Equivalent shear strain rate [1/s]
 γ_{ij} Shear strain rate tensor [1/s]
 ρ Specific mass [kg/m^{3}]
 τ Local shear stress magnitude [Pa]
 τ^{*} Dimensionless local shear stress: τ^{*} = τ/τ_{0}
 τ_{0} Yield stress [Pa]
 τ_{ij} Shear (viscous) stress tensor [Pa]
 Γ_{i} Perimeter of the physical domain [i = 1, 2, 3, 4]
 Ω_{AUR} Dimensionless unyielded area: Ω_{AUR} = A_{AUR}/H^{2}
Publication Dates

Publication in this collection
09 Dec 2019 
Date of issue
JulSep 2019
History

Received
01 Aug 2018 
Reviewed
13 Dec 2018 
Accepted
29 Mar 2019