NUMERICAL STUDY OF VISCOPLASTIC FLOW IN A T-BIFURCATION: IDENTIFICATION OF STAGNANT REGIONS

- 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 T-bifurcation channel is adopted to illustrate the procedure. The rheological behavior of Carbopol ® 980 was simulated using the Herschel-Bulkley 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.


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, 2017). 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., 1983), so that both yielded and unyielded regions exist in the flow. The ideal viscoplastic models, such as Bingham plastic, Herschel-Bulkey and Casson, are discontinuous and numerical solutions for complex geometries require regularized models, such as the biviscosity equation (Tanner and Milthorpe, 1983) and models proposed by Bercovier and Engelman (1980) and Papanastasiou (1987). The latter was adopted in this study to investigate the hydrodynamic behavior of a Herschel-Bulkley fluid in a T-bifurcation. The Papanastasiou regularization uses an exponential modification now widely applied in numerical studies involving viscoplastic materials (Mitsoulis, 2007).
In unyielded regions, the stress state falls below the yield stress threshold and the viscoplastic material behaves as a very high-viscosity 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., 1992;Mitsoulis et al., 1993). Papanastasiou and Boudouvis (1997) evaluated (a) the non-shear 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., 2007). 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, 2007). Carbopol ® dispersions are found in many everyday products, from toothpastes to tile cleaners, and are also useful vehicles for functional ingredients (Roberts and Barnes, 2001). As a thickening solution used in the cosmetic industry, it has been largely applied in preparation of gels (Corrêa, 2005), with the purpose of emulsification, stabilization and rheological control (Kim, 2003).
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 T-bifurcation 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 T-bifurcation. 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, 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., 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., 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 Herschel-Bulkley model (Rudert and Schwarze, 2009). The literature also shows that thixotropic characteristics of these gels can, in principle, be neglected (Coussot, 2014). 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 Herschel-Bulkley) do not present limitations to analytical solutions of simple problems (Bird et al., 1983), 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.  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., 1997). Papanastasiou (1987) 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: 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 cell-centered formulation and a second order scheme for discretization of the viscous terms. The SIMPLE method was adopted to approach the pressure-velocity coupling. The odd-even decoupling was resolved using the Rhie-Chow interpolation scheme. The viscosity function was determined using the regularized Herschel-Bulkey-Papanastasiou 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 two-dimensional plane channel using ANSYS FLUENT ® (based on the finite volume method) was compared against solutions (i) obtained using an in-house code developed by the authors based on finite differences (FDM) and (ii) reported by Boualit et al. (2011) 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 Bingham-Papanastasiou 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); non-slip condition at the channel walls (Γ 2 : u * = 0, v * = 0); fully-developed 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 ).
In order to guarantee the fully-developed velocity distribution towards the exit, a channel length L = 80H where η 0 is a constant viscosity parameter and τ 0 is the yield stress. The Papanastasiou regularization applied to the Herschel-Bulkley 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 was considered. A mesh containing 300 x 30 elements was also used in the simulations. A comparative assessment of the axial velocities for fully-developed flow for different Bingham numbers, Bn, against the results reported by Boualit et al (2011) 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. (2011) and the analytical study carried out by Quaresma and Macêdo (1998). 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. (2011).
Flow in a T-bifurcation channel -mesh refinement study: The distribution process of a viscoplastic material is simulated considering a T-bifurcation 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 value of the regularization parameter m was defined according to the works of , , Alexandrou et al. (2001) and Boualit et al. (2011). The Reynolds number for the Herschel-Bulkley viscoplastic fluid can be derived as:

Analysis of the flow in a channel with T-bifurcation:
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 T-bifurcation 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  Assessment of mesh dependency was evaluated by the dimensionless pressure drop, Dp * , and fullydeveloped dimensionless velocities u * (x = 15H) and v * (y = ±20H) 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 = ±20H) 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 fully-developed regions, the mesh containing 28,125 elements was selected for the study.
the work of Zisis and Mitsoulis (2002). Therefore, by using Equations (5) and (6), the magnitude of the local shear stress is so that, if t < t 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, t * = t/t 0 , is clipped at t * < 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 T-bifurcation. 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.   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 t * < 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 Herschel-Bulkley-Papanastasiou 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. (2001), who also found similar results by evaluating a Herschel-Bulkley 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 Herschel-Bulkley fluids. The (12) 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, W (i) AUR , where t * < 1, as: 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 T-type bifurcation were evaluated by numerical solution. The finite volume method was applied using the ANSYS FLUENT ® package associated with Papanastasiou's regularization for Herschel-Bulkley 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 T-bifurcation channel addressed in this work by providing the smallest fraction of all types of unyielded stagnant regions. 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, W 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 T-bifurcation channel. This