The Role of Elasticity in the Vortex Formation in Polymeric Flow around a Sharp Bend

Featured Application: Elasticity is predicted to play a signiﬁcant role in the formation of vortices in polymeric ﬂow in a sharp bend. The polymer dilution and the ﬂow rate determine if the formation of the vortices occurs upstream or downstream of the bend corner. Abstract: Fluid dynamic simulations using the FENE-P model of polymer physics are compared to those of an incompressible Newtonian ﬂuid base case in order to understand the role of elasticity in the formation of vortices in a 90 ◦ bend narrow channel. The analysis bridges the ﬂow behavior of a purely elastic ﬂuid and that of a Newtonian ﬂuid. We evaluated how four dimensionless numbers—Reynolds number (Re), Weissenberg number (Wi), viscosity ratio ( β ), and elasticity number (El)—affect the formation of vortices. It is shown that increasing Re and Wi, or lowering β will cause vortices to grow in size. Two phase space diagrams, β vs. El and β vs. Re, were created to show the range of values where inertial and elastic vortices form. Both diagrams have three zones. Depending on the polymer viscosity ratio and the elasticity number, the vortices form either upstream of the bend (elasticity driven) or form downstream of the bend (inertia driven), are suppressed. Our predictions are in good agreement with previous experimental and numerical works. vortices for polymer solutions with speciﬁc rheological properties and ﬂow conditions. Our


Introduction
Many modern-day products, such as rubber tires and plastic bags, are made from polymers. These items are made by specialized processing machines that handle polymer solutions in pipes, conduits, and accessories such as bends and elbows. These polymeric solutions, as with any other liquid, are subjected to many fluid dynamic effects. One such effect is the formation of recirculation zones, or vortices in abrupt, changes of flow direction as in contractions and bends. These vortices are a feature where some amount of the fluid becomes trapped in a cyclone like structure near or around corners, justifying those authors that refer to them as separation bubbles [1]. Recirculation zones can have detrimental effects on the flow of polymers and polymer solutions, affecting their manufacturing processes. As a result, it is important to know the conditions under which vortices form. The presence of vortices near or around corners is one of the most significant alterations in channel flow for both non-Newtonian and Newtonian fluids. Vortices are well known to form in Newtonian fluids due to the effect of inertia [2]. However, for polymers and polymer solutions not all the physics involved are fully understood [3,4].
Understanding the flow behavior of polymers and polymer solutions is essential to design and optimize fluid flow systems in many practical applications. For example, in microfluidics the polymer concentration has a significant effect on viscoelastic behavior by altering the base flow or result in flow instabilities. In that regard, Gulati et al. [5] studied the flows of dilute and semi-dilute polymer solutions in sharp 90 • micro-bends in channels of rectangular cross-section. Their flow visualizations show that a vortex is present in the inner, upstream corner of the bend and grows with increasing Reynolds and Weissenberg numbers for flows of shear-thinning, semi-dilute polymeric solutions. They reported that secondary flows were not present for Newtonian flows under similar conditions and that a vortex is absent for flow of a dilute, non-shear thinning PEO solution. They concluded that shear-thinning appears to be central to the presence of an elastic secondary flow in this geometry. Their experiments were carried at very low Reynolds number (10 −6 < Re < 0.03), and Weissenberg numbers ranging from 0.42 to 126.
More recently, Kim et al. [6] reported instabilities in viscoelastic flow in a 90 • bent channel. They observed that the flow instability in an aqueous PEO solution occurs when the concentration of PEO is as low as 50 ppm. Investigating the effects of the polymer concentration, flow rate, and elasticity number, they found that the flow is stabilized in shear-thinning fluids, whereas the flow instability is amplified when both elastic and inertial effects are pronounced. Their experiments were carried out at Reynolds numbers ranging from 0.3 to 3.0 and Weissenberg numbers of 0 (Newtonian flow) to 40.
Both studies [5,6] coincide in pointing to the shear thinning properties of the solution, a decrease in viscosity under shear strain, as the main reason for the formation of vortices and secondary flows.
At the macroscopic level, it is also well known that polymers, as well as surfactants, are frequently added to Newtonian fluids with the purpose of reducing friction losses in straight pipe turbulent flow. In some cases, the drag-reducing rate is as high as 75% [7] making them attractive for using in complex industrial pipe flow systems. Friction losses in industrial piping systems are mostly due to accessories such as bends, tees, and valves, rather than in the straight pipes. Understanding how polymers affect the Newtonian flow in bends and accessories would help in evaluating their drag reduction potential and application in intricate piping systems.
Munekata et al. [8] studied the bend flow characteristics of two surfactant solutions, experimentally and numerically. They found that the drag, friction coefficient, increases or reduces depending on the average bulk velocity of the solution. As the solution concentration increases a larger bulk velocity is required to observe a reduction in the drag. Although the drag reduction's effects of the surfactants in the solution are lower in bend flow than in a straight pipe, the authors attribute the drag reduction to the suppression of the centrifugal effect and reducing the secondary flow due to the viscoelastic properties of the solutions (normal stress effect). Their analysis was performed at a high Reynolds number for both a Newtonian fluid and a highly viscoelastic fluid. They observed and predicted smaller velocity gradients near the wall for viscoelastic fluid flow than that for Newtonian fluid flow.
Even in purely Newtonian fluid flow, understanding the hydrodynamic behavior in a bend channel or a pipe elbow remains relevant today. Matsumoto et al. [1] investigated the flow dynamics for a Newtonian fluid in a bent channel via two-dimensional direct numerical simulations. They investigated the flow structure along the channel as a function of both the bend angle and the Reynolds number. Their numerical work suggests a scaling relation between the shape of the separation bubble, a downstream vortex after the bend corner, and the flow conductance. Their simulations were carried at high Reynolds number but only for a Newtonian fluid. Nevertheless, they present an integrated phase diagram for the flow dynamics, where depending on bend angle and Reynolds number either the flow is uniform with no recirculation vortex forming, a vortex forms downstream of then bend corner, or vortices shed intermittently from the bend corner.
In this paper we address the role of elasticity in the formation of vortices. With our analysis we are bridging the flow behavior of a polymeric fluid and that of a Newtonian fluid in a 90 • bend narrow channel.
In agreement with previous works [5,6], we predict the formation of elastic and inertial vortices for polymer solutions with specific rheological properties and flow conditions. Our primary interest is exploring the possibility of controlling the alteration of the flow, either suppressing or promoting a vortex or separation bubble, by modifying the underlying properties of the polymer solution.
The remainder of this paper is organized as follows. In Section 2, we introduce our hydrodynamic and polymer models. In Section 3, we describe in detail our numerical procedure including simulation parameters, and the computational domain including a few details about the meshing. In Section 4, we present the numerical results, specifically comparing Newtonian and polymer flow and discussing the role of elasticity and Reynolds number on vortex formation. We also construct two flow phase diagrams. In Section 5, we summarize our results and briefly comment on the future perspectives of the present work.

Governing Equations
Two types of fluids were studied in this research, a Newtonian fluid and a polymeric fluid. Newtonian flows are characterized by Navier-Stokes equations. Polymer flows require additional equations to characterize the elastic behavior of the polymer chains. In this work, the polymer chains are characterized using the FENE-P model (finitely extensible nonlinear elastic model with Peterlin closure [9]). The FENE-P model was chosen due to its versatility and simplicity in characterizing polymer behavior. In the FENE-P model, polymers are represented as dumbbells, two masses connected by a spring, as in many other polymer models [10]. Nevertheless, in the FENE-P model the spring has a finite stretching limit.
For this work it is assumed that the fluid is incompressible for both Newtonian and polymer-based flows, for which the continuity equation becomes: where v is the velocity of the fluid. The general momentum equation of motion is given by: where Π is the total stress tensor: with τ as the extra stress tensor, p the hydrostatic pressure, and I the identity tensor.
For an incompressible Newtonian fluid, the stress tensor τ is given by: .
where η s is the Newtonian shear viscosity, and . γ is the strain rate tensor. Combining (3)-(5), substituting into (2), and neglecting the effects of gravity yields the well-known Navier-Stokes equation of fluid dynamics, see below. For a detailed derivation of (6) see ref. [11].
For complex fluids such as polymers and polymer solutions, Equation (4) is not sufficient to describe the dynamics of the flows. An additional stress term is added to the Newtonian constitutive equation. For complex fluids: Appl. Sci. 2021, 11, 6588 4 of 19 where τ N denotes the Newtonian stress component, and τ P denotes the stress component characterizing the polymer behavior. The equation used to estimate τ P depends on the particular model of polymer physics being employed. For the purposes of this research, the FENE-P constitutive equation is given by [10]: where G O is the elastic modulus, and λ is the relaxation time of the polymers. The conformation tensor A measures the stretch and orientation of the polymers. The upper-convective derivative [12], denoted by ( ) (1) , is given by: and: with b representing the square of the maximum extensibility. The total stress in the polymer solution is then obtained by combining Equations (7) and (8), and substituting that result back into (3):

Numerical Solver
The numerical solver implemented to solve the governing equations for the present work is the open-source programming package RheoTool. RheoTool is an open-source toolbox based on the OpenFOAM ® library to simulate the flow of Generalized Newtonian Fluids (GNF) and viscoelastic fluids [13].
RheoTool is a modification of the viscoelastic solver available in the OpenFOAM ® toolbox [14]. The main goal of the modification was to improve its stability for differential-type constitutive equations. The major contributions of RheoTool are using the log-conformation approach to solve Oldroyd-B type constitutive equations, handling high-resolution schemes with a componentwise and deferred correction approach to discretize the convective terms, and introducing a new stress-velocity coupling term together with the well-known SIM-PLEC algorithm for pressure-velocity coupling [15].
In the present work, the polymer is characterized by the FENE-P model which is solved using the log-conformation approach [16,17], and selecting the stress-velocity coupling as the stabilization method.
In terms of discretization, gradient terms are discretized using the Gauss scheme with linear interpolation, Laplacian terms are discretized using the Gauss scheme (only choice) with linear interpolation and the corrected scheme for the surface normal gradient, and the convective terms are discretized using the CUBISTA scheme [18].
The solver chosen for all (asymmetric) equations is the Preconditioned (bi-) Conjugate Gradient with the Diagonal incomplete-Cholesky (LU) preconditioner.
The solution is advanced in time using the Euler scheme with adjustable time steps. A minimum time step of 10 −6 s is set together with a maximum Courant number of 2.0 and a maximum time step of 10 −2 s. This time step condition allows for speed up of the convergence of transient simulations to steady state.

Simulation Parameters
The main goal of this research is to evaluate the role that elasticity plays in the formation of vortices in flow around a sharp corner. For that purpose, we performed Appl. Sci. 2021, 11, 6588 5 of 19 multiple simulations to create two phase diagrams in order to show the overall effects of four dimensionless parameters, Wi, Re, El, and β. These parameters are defined as follows:

•
The Weissenberg number Wi = λ·U/H, where λ is the relaxation time of the polymer, U is the characteristic velocity, and H is the characteristic length scale of the geometry (channel height). The Weissenberg number is a ratio of the polymeric timescale to a convective timescale, akin to an elastic to viscous forces ratio in the context of this work.

•
The Reynolds number Re = ρ·U·H/η 0 , where ρ and η 0 are, respectively, the density and the zero-shear rate viscosity of the solution. The Reynolds number is a measure of the ratio of inertia to viscous forces.

•
The elasticity number, a derived parameter, characterizes the balance of elastic and inertial forces in the fluid and is defined as El = Wi/Re.

•
The solvent viscosity ratio β = η s /η o , where η s is the Newtonian solvent viscosity and η o is the zero-shear rate viscosity of the solution.
An important note about the elasticity number El and the viscosity ratio β is that they are both material parameters that depend on η o , and are independent of the fluid velocity. The zero-shear rate viscosity, η o , physically corresponds to the concentration of polymers in the solution. A decreasing El indicates a decrease in the polymer concentration. Conversely, an increase in β corresponds also to a decrease in the polymer concentration, which is the parameter typically varied during experimental work.
Two sets of simulations were conducted in order to meet this goal. The first group of simulations conducted during this research were used to create a Weissenberg number versus Reynolds number phase space diagram. The Weissenberg numbers evaluated were 0.5, 1, and 1.5. Simulations were run with viscosity ratios, β values, ranging from 0 to 1 in 0.1 increments. A β value of 0 represents a polymer melt. A β value of 1 represents the viscosity of a Newtonian fluid (solvent). A second group of simulations was conducted in order to create a viscosity ratio, β, versus elasticity number, El, phase diagram. All of these simulations were performed with a constant value for the Weissenberg number (Wi = 1). Figure 1 shows a representation of the computational domain used to model the right-angle bend geometry. The fluid flows into the system at the top left, where a uniform inlet velocity (for a given flowrate Q) is imposed, travels around the sharp bend, and then exits out at the bottom, where the outlet pressure is set. At the walls, no-slip and no-penetration conditions are set for the velocities, and zero gradient for the pressure and the stress tensors. The bend used in the simulations has a ten to one ratio of length to channel height, H, in both branches. This ratio is more than enough to obtain a fully developed flow well upstream of the corner bend and for the exit, being far enough away to prevent any upstream effect from the outlet boundary conditions, see references [19,20], and also Appendix A. The geometry used follows the experimental work of Gulati et al. [5] and is typical of benchmark cases and similar configurations that have been evaluated both numerically and experimentally, see [1,6,21,22]. The bend used in the simulations has a ten to one ratio of length to channel height, H, in both branches. This ratio is more than enough to obtain a fully developed flow well upstream of the corner bend and for the exit, being far enough away to prevent any upstream effect from the outlet boundary conditions, see references [19,20], and also Appendix A. The geometry used follows the experimental work of Gulati et al. [5] and is typical of benchmark cases and similar configurations that have been evaluated both numerically and experimentally, see [1,6,21,22]. Figure 2 presents an enlarged portion of the bend corner geometry with the progression in mesh refinement used for the simulations. Preliminary runs were performed in the coarse mesh shown in Figure 2a. The mesh was halved twice for the final run of every simulation, Figure 2c. An even finer fourth mesh was used in a few simulations to estimate grid convergence regarding the size predicted for inertial and elastic vortices including the phase maps.

Meshing
The bend used in the simulations has a ten to one ratio of length to channel height, H, in both branches. This ratio is more than enough to obtain a fully developed flow well upstream of the corner bend and for the exit, being far enough away to prevent any upstream effect from the outlet boundary conditions, see references [19,20], and also Appendix A. The geometry used follows the experimental work of Gulati et al. [5] and is typical of benchmark cases and similar configurations that have been evaluated both numerically and experimentally, see [1,6,21,22]. Figure 2 presents an enlarged portion of the bend corner geometry with the progression in mesh refinement used for the simulations. Preliminary runs were performed in the coarse mesh shown in Figure 2a. The mesh was halved twice for the final run of every simulation, Figure 2c. An even finer fourth mesh was used in a few simulations to estimate grid convergence regarding the size predicted for inertial and elastic vortices including the phase maps.      Figure 3b compares the spanwise velocity component for the same four grids. The weighted average deviation error outside the vortex is 14% and 26% inside the vortex, respectively (in both cases, note the very low actual value of the predicted velocity). Lastly, Figure 3c compares the axial stress. The weighted average devia-  inside the vortex. Similarly, Figure 3b compares the spanwise velocity component for the same four grids. The weighted average deviation error outside the vortex is 14% and 26% inside the vortex, respectively (in both cases, note the very low actual value of the predicted velocity). Lastly, Figure 3c compares the axial stress. The weighted average deviation error outside the vortex is 4.0% and 8.3% inside the vortex, respectively.

Meshing
Similar analyses were conducted for many cases, including those with inertial vortices located downstream of the bend corner. In a typical high Reynolds number case, Re = 100, the weighted average deviation errors of all variables are significantly smaller. For the streamwise velocity component, the weighted average deviation error outside the vortex for the fine grid is 0.02% and 0.4% inside the vortex, respectively. For the spanwise velocity component, the weighted average deviation error outside the vortex is 0.09% and 0.5% inside the vortex, respectively. Additionally, for the strain rate, the weighted average deviation error outside the vortex is 0.04% and 0.28% inside the vortex, respectively.
In general terms, it was found that predicting the location and length of the vortex is the most significant challenge in attaining grid convergence. Outside the vortex and far from it, the grid convergence index (GCI) for the fine grid is excellent, ranging from 0.0 to 1.5% in all variables. Within the vortex, the GCI for the fine grid is acceptable, ranging from 0.1% to 25% in all variables. Again, it should be noted that actual values for the spanwise velocity and the normal stress within the vortex are significantly smaller when compared to the streamwise velocity and the axial stresses, by two to three orders of magnitude. The mesh and grid convergence analyses were conducted following the Procedure for Estimation and Reporting of Uncertainty Due to Discretization in CFD Applications by Roache et al. [23].
The formation of vortices in polymer flow has been studied extensively in planar contractions and the size and strength have been shown to be strongly dependent on the mesh [15,24,25]. For example, the effect of mesh refinement on the size of the vortices was highlighted by Alves et al. [26] when developing benchmark solutions for the flow of Oldroyd-B and PTT fluids in planar contractions. They found that a very high degree of mesh fineness was required to obtain accurate results with the Oldroyd-B fluid, while the PTT fluid in general did not require the finest meshes. These authors relied on Richardson's extrapolation to measure the level of convergence.

Newtonian Flow
To better understand the effects the elasticity of polymers have on the formation of vortices in a right angle bend geometry, simulations of Newtonian fluid flow were conducted first in order to obtain a baseline or reference. It is well known that in Newtonian flow vortices form downstream of a bend corner at sufficiently high Reynolds number [22,27].
To establish the baseline, the Reynolds number value at which the vortices first appear was sought. This value, called Re crit for the duration of this paper, was found to be approximately 30. For Reynolds number values larger than 30, a clear vortex forms downstream of the bend corner. This critical Reynolds number value is in line with other numerical work [1,2].
For Newtonian fluids, the Reynolds number is the single parameter required to describe the flow. It is expected that higher Reynolds numbers should lead to larger vortices downstream of the bend corner. We refer to these vortices forming downstream as inertial vortices from now on. Figure 4 presents two inertial vortices. The vortex in Figure 4a  For Newtonian fluids, the Reynolds number is the single parameter required to describe the flow. It is expected that higher Reynolds numbers should lead to larger vortices downstream of the bend corner. We refer to these vortices forming downstream as inertial vortices from now on. Figure 4 presents two inertial vortices. The vortex in Figure 4a corresponds to a flow with a Reynolds number of 50, while Figure 4b corresponds to Re = 100.  Comparing both figures, it can be seen that the second vortex at Re = 100 is predicted to be significantly larger than the vortex that forms at Re = 50. This comparison shows that a higher Reynolds number value leads to a larger inertial vortex in Newtonian flow. It should be noted that the size (length) of the vortices in these and all figures that follow was determined numerically by the change in direction of the tangential velocity near the wall, equivalent to a change in sign of the shear stress (skin friction), see Appendix B.
To validate both, the statement that the inertial vortex size increases with Reynolds number and that 'RheoTool' is a valid numerical open source tool for the solution of both Newtonian and non-Newtonian flow, we modeled the bend Newtonian flow using the commercial tool ANSYS Fluent v19.1. Figure 5 compares the non-dimensional vortex length (vortex length to channel height ratio) predicted with both Fluent and Rheotool up to a Reynolds number of 250. The values predicted by RheoTool are within 1% of those predicted by Fluent. The figure clearly shows that the vortex length (l) increases with Reynolds number. Additionally, the comparison also indicates that RheoTool is as good as Fluent in solving the Navier-Stokes equations for Newtonian flow. It should be noted Comparing both figures, it can be seen that the second vortex at Re = 100 is predicted to be significantly larger than the vortex that forms at Re = 50. This comparison shows that a higher Reynolds number value leads to a larger inertial vortex in Newtonian flow. It should be noted that the size (length) of the vortices in these and all figures that follow was determined numerically by the change in direction of the tangential velocity near the wall, equivalent to a change in sign of the shear stress (skin friction), see Appendix B.
To validate both, the statement that the inertial vortex size increases with Reynolds number and that 'RheoTool' is a valid numerical open source tool for the solution of both Newtonian and non-Newtonian flow, we modeled the bend Newtonian flow using the commercial tool ANSYS Fluent v19.1. Figure 5 compares the non-dimensional vortex length (vortex length to channel height ratio) predicted with both Fluent and Rheotool up to a Reynolds number of 250. The values predicted by RheoTool are within 1% of those predicted by Fluent. The figure clearly shows that the vortex length (l) increases with Reynolds number. Additionally, the comparison also indicates that RheoTool is as good as Fluent in solving the Navier-Stokes equations for Newtonian flow. It should be noted that this Newtonian inertial vortex would become unstable at a sufficiently high Reynolds number value [1].
Appl. Sci. 2021, 11, x FOR PEER REVIEW that this Newtonian inertial vortex would become unstable at a sufficiently high R number value [1].

Polymer Flow
As previously discussed, there are four dimensionless numbers-Wi, Re, El, typically used to characterize polymer fluid dynamics. Nevertheless, the elastici ber (El) is defined as the ratio of the Weissenberg and Reynolds numbers. Conse only three dimensionless parameters are effectively needed to characterize polym Evaluating and predicting the effect on these parameters on polymer flow and vo mation is the main goal of this research.
It has already been shown, for the Newtonian case, that a difference in R number means a difference in vortex size. For polymer fluids, even though elas present, this should be no different. Holding the Weissenberg number constant c ing two simulations with different elasticity numbers should give insight into h Reynolds number affects the vortices.

Inertial Vortices
The vortices in Figure 6 correspond to flows at different Reynolds numbers same Weissenberg number of 1.5 and a viscosity ratio of 0.1. Larger downstream are predicted as the Reynolds number is increased. With both Wi and β held con should be this variation in Reynolds number that is accountable for the increase i size. It should be noted that increasing the Reynolds number while holding both β constant is equivalent to reducing the elasticity number.

Polymer Flow
As previously discussed, there are four dimensionless numbers-Wi, Re, El, and βtypically used to characterize polymer fluid dynamics. Nevertheless, the elasticity number (El) is defined as the ratio of the Weissenberg and Reynolds numbers. Consequently, only three dimensionless parameters are effectively needed to characterize polymer flow. Evaluating and predicting the effect on these parameters on polymer flow and vortex formation is the main goal of this research.
It has already been shown, for the Newtonian case, that a difference in Reynolds number means a difference in vortex size. For polymer fluids, even though elasticity is present, this should be no different. Holding the Weissenberg number constant comparing two simulations with different elasticity numbers should give insight into how the Reynolds number affects the vortices.

Inertial Vortices
The vortices in Figure 6 correspond to flows at different Reynolds numbers for the same Weissenberg number of 1.5 and a viscosity ratio of 0.1. Larger downstream vortices are predicted as the Reynolds number is increased. With both Wi and β held constant, it should be this variation in Reynolds number that is accountable for the increase in vortex size. It should be noted that increasing the Reynolds number while holding both Wi and β constant is equivalent to reducing the elasticity number.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 9 of 20 that this Newtonian inertial vortex would become unstable at a sufficiently high Reynolds number value [1].

Polymer Flow
As previously discussed, there are four dimensionless numbers-Wi, Re, El, and βtypically used to characterize polymer fluid dynamics. Nevertheless, the elasticity number (El) is defined as the ratio of the Weissenberg and Reynolds numbers. Consequently, only three dimensionless parameters are effectively needed to characterize polymer flow. Evaluating and predicting the effect on these parameters on polymer flow and vortex formation is the main goal of this research.
It has already been shown, for the Newtonian case, that a difference in Reynolds number means a difference in vortex size. For polymer fluids, even though elasticity is present, this should be no different. Holding the Weissenberg number constant comparing two simulations with different elasticity numbers should give insight into how the Reynolds number affects the vortices.

Inertial Vortices
The vortices in Figure 6 correspond to flows at different Reynolds numbers for the same Weissenberg number of 1.5 and a viscosity ratio of 0.1. Larger downstream vortices are predicted as the Reynolds number is increased. With both Wi and β held constant, it should be this variation in Reynolds number that is accountable for the increase in vortex size. It should be noted that increasing the Reynolds number while holding both Wi and β constant is equivalent to reducing the elasticity number. It can then be stated that larger Reynolds numbers lead to larger inertial vortices in polymer flow. Conversely, reducing the Reynolds number by increasing the elasticity of the polymer solution should reduce the inertial vortex size.
Like Newtonian flow, current predictions indicate that larger Reynolds numbers in polymer flow should lead to larger inertial vortex sizes. This is consistent with laminar Newtonian flow physics, but what role do the polymers play in inertial vortex sizes? A way to evaluate the effect of the polymers on the flow would be to compare a Newtonian inertial vortex to a polymer vortex at the same Reynolds number.  Figure  7d corresponds to a polymer flow with Wi = 1.5 (El = 0.015). All three polymer inertial vortices are predicted to be comparable in size, within 5%; with the one corresponding to a Weissenberg number of 1.0 being the only one significantly larger than its Newtonian counterpart, but only about 10%, just slightly above the average grid convergence error of 8%. This is an indication that adding polymers to the flow at this high Reynolds number is predicted to have a slight to nonexistent effect on the fluid dynamics of the flow around a bend. The Weissenberg number effect on inertial vortices at high Reynolds numbers is predicted to only be moderate. It can then be stated that larger Reynolds numbers lead to larger inertial vortices in polymer flow. Conversely, reducing the Reynolds number by increasing the elasticity of the polymer solution should reduce the inertial vortex size.
Like Newtonian flow, current predictions indicate that larger Reynolds numbers in polymer flow should lead to larger inertial vortex sizes. This is consistent with laminar Newtonian flow physics, but what role do the polymers play in inertial vortex sizes? A way to evaluate the effect of the polymers on the flow would be to compare a Newtonian inertial vortex to a polymer vortex at the same Reynolds number.  Figure 7d corresponds to a polymer flow with Wi = 1.5 (El = 0.015). All three polymer inertial vortices are predicted to be comparable in size, within 5%; with the one corresponding to a Weissenberg number of 1.0 being the only one significantly larger than its Newtonian counterpart, but only about 10%, just slightly above the average grid convergence error of 8%. This is an indication that adding polymers to the flow at this high Reynolds number is predicted to have a slight to nonexistent effect on the fluid dynamics of the flow around a bend. The Weissenberg number effect on inertial vortices at high Reynolds numbers is predicted to only be moderate. It can then be stated that larger Reynolds numbers lead to larger inertial vortices in polymer flow. Conversely, reducing the Reynolds number by increasing the elasticity of the polymer solution should reduce the inertial vortex size.
Like Newtonian flow, current predictions indicate that larger Reynolds numbers in polymer flow should lead to larger inertial vortex sizes. This is consistent with laminar Newtonian flow physics, but what role do the polymers play in inertial vortex sizes? A way to evaluate the effect of the polymers on the flow would be to compare a Newtonian inertial vortex to a polymer vortex at the same Reynolds number.  Figure  7d corresponds to a polymer flow with Wi = 1.5 (El = 0.015). All three polymer inertial vortices are predicted to be comparable in size, within 5%; with the one corresponding to a Weissenberg number of 1.0 being the only one significantly larger than its Newtonian counterpart, but only about 10%, just slightly above the average grid convergence error of 8%. This is an indication that adding polymers to the flow at this high Reynolds number is predicted to have a slight to nonexistent effect on the fluid dynamics of the flow around a bend. The Weissenberg number effect on inertial vortices at high Reynolds numbers is predicted to only be moderate. The contour plots in Figure 8 show two inertial vortices both corresponding to a moderate Reynolds number of 50 also with β = 0.1. The first vortex, Figure 8a, has a Weissenberg number of 0.5, while the second, Figure 8b, corresponds to Wi = 1.0. Again, these inertial polymer vortices are comparable in size, in this case within 2%, indicating that the influence of the polymer in the size of the inertial vortex is also predicted to be moderate at a moderate Reynolds number.   The contour plots in Figure 8 show two inertial vortices both corresponding to a moderate Reynolds number of 50 also with β = 0.1. The first vortex, Figure 8a, has a Weissenberg number of 0.5, while the second, Figure 8b, corresponds to Wi = 1.0. Again, these inertial polymer vortices are comparable in size, in this case within 2%, indicating that the influence of the polymer in the size of the inertial vortex is also predicted to be moderate at a moderate Reynolds number. The contour plots in Figure 8 show two inertial vortices both corresponding to a moderate Reynolds number of 50 also with β = 0.1. The first vortex, Figure 8a, has a Weissenberg number of 0.5, while the second, Figure 8b, corresponds to Wi = 1.0. Again, these inertial polymer vortices are comparable in size, in this case within 2%, indicating that the influence of the polymer in the size of the inertial vortex is also predicted to be moderate at a moderate Reynolds number.     From Figure 9 it could be inferred that adding polymers to the solution does not affect the size of the inertial vortices. Figure 10 presents the same data as Figure 9 but as a function of the elasticity number. At every Weissenberg number, the predictions clearly indicate that increasing elasticity leads to a rapid decrease in vortex size. This figure also highlights that it is the Reynolds number, indeed the inertial effects, that is predicted to drive the size of the inertial vortex, with the Weissenberg number playing a secondary role. This figure also seems to indicate that there is a maximum size for the inertial vortex at any given Reynolds number. This maximum vortex size is predicted for the cases where the inertial and elastic forces are balanced: Weissenberg equal to 1.0.

Elastic Vortices
As the Reynolds number is reduced, the downstream or inertial vortices are predicted to disappear. However, as the Reynolds number is reduced further, vortices are predicted to form upstream of the corner bend. We refer to these vortices forming upstream of the bend corner as elastic vortices. At low Reynolds numbers, it is the Weissenberg number that is predicted to determine the size of the elastic vortex. To evaluate how changing the Weissenberg number affects vortices that form upstream of the bend corner, Figure 11 compares an elastic vortex for Wi = 1 and El = 1·10 6 to a vortex for Wi = 1.5 and From Figure 9 it could be inferred that adding polymers to the solution does not affect the size of the inertial vortices. Figure 10 presents the same data as Figure 9 but as a function of the elasticity number. At every Weissenberg number, the predictions clearly indicate that increasing elasticity leads to a rapid decrease in vortex size. This figure also highlights that it is the Reynolds number, indeed the inertial effects, that is predicted to drive the size of the inertial vortex, with the Weissenberg number playing a secondary role. This figure also seems to indicate that there is a maximum size for the inertial vortex at any given Reynolds number. This maximum vortex size is predicted for the cases where the inertial and elastic forces are balanced: Weissenberg equal to 1.0. From Figure 9 it could be inferred that adding polymers to the solution does not affect the size of the inertial vortices. Figure 10 presents the same data as Figure 9 but as a function of the elasticity number. At every Weissenberg number, the predictions clearly indicate that increasing elasticity leads to a rapid decrease in vortex size. This figure also highlights that it is the Reynolds number, indeed the inertial effects, that is predicted to drive the size of the inertial vortex, with the Weissenberg number playing a secondary role. This figure also seems to indicate that there is a maximum size for the inertial vortex at any given Reynolds number. This maximum vortex size is predicted for the cases where the inertial and elastic forces are balanced: Weissenberg equal to 1.0.

Elastic Vortices
As the Reynolds number is reduced, the downstream or inertial vortices are predicted to disappear. However, as the Reynolds number is reduced further, vortices are predicted to form upstream of the corner bend. We refer to these vortices forming upstream of the bend corner as elastic vortices. At low Reynolds numbers, it is the Weissenberg number that is predicted to determine the size of the elastic vortex. To evaluate how changing the Weissenberg number affects vortices that form upstream of the bend corner, Figure 11 compares an elastic vortex for Wi = 1 and El = 1·10 6 to a vortex for Wi = 1.5 and

Elastic Vortices
As the Reynolds number is reduced, the downstream or inertial vortices are predicted to disappear. However, as the Reynolds number is reduced further, vortices are predicted to form upstream of the corner bend. We refer to these vortices forming upstream of the bend corner as elastic vortices. At low Reynolds numbers, it is the Weissenberg number that is predicted to determine the size of the elastic vortex. To evaluate how changing the Weissenberg number affects vortices that form upstream of the bend corner, Figure 11 compares an elastic vortex for Wi = 1 and El = 1·10 6 to a vortex for Wi = 1.5 and El = 1.5·10 6 . Both simulations correspond to a viscosity ratio of 0.1 and a Reynolds number of 1·10 −6 . It is clear that the two vortices are predicted to be of different sizes, with the vortex for Wi = 1.5 being larger. Since the elasticity number and the Weissenberg number are linked, the vortices are predicted to be different in size, most likely due to the difference in Weissenberg numbers.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 13 of 20 El = 1.5·10 6 . Both simulations correspond to a viscosity ratio of 0.1 and a Reynolds number of 1·10 −6 . It is clear that the two vortices are predicted to be of different sizes, with the vortex for Wi = 1.5 being larger. Since the elasticity number and the Weissenberg number are linked, the vortices are predicted to be different in size, most likely due to the difference in Weissenberg numbers.
(a) (b) Lastly, we evaluate the effect of the viscosity ratio. As explained previously, this ratio can take any value between 0 and 1. Since a β value of 1 defines a purely Newtonian fluid, we expect that fluids with viscosity ratios less than 1 should exhibit polymer behavior similar to increasing the Weissenberg number in polymer flow. The elastic vortex should therefore be larger for smaller β values because of its reliance on elasticity to form.
The elastic vortices in Figure 12 show the effect of decreasing β; from β = 0.4 to β = 0.1. The second elastic vortex, lower viscosity ratio, is larger than expected, indeed, significantly larger, roughly three times. Reducing β is predicted to have a similar effect to that of increasing the Weissenberg number, increasing the size of the elastic vortex. However, the presence of elasticity is not sufficient to produce elastic vortices; it was found that elastic vortices are only predicted to form for viscosity ratios not exceeding a certain value. In other words, a dilute polymer solution, mostly solvent with a viscosity ratio near 1, is always predicted to have Newtonian-like flow behavior without a vortex forming upstream of the bend corner. Conversely, a concentrated polymer solution, mostly polymer with a viscosity ratio near 0, is predicted to have secondary flows, vortices upstream of the bend corner, for a large range of Reynolds numbers. These predictions are in line with the experiments of Gulati [5] who found that stable elastic vortices form within the flow of semidilute DNA solutions in a 90° micro bend channel above a certain threshold for the solution elasticity. Lastly, we evaluate the effect of the viscosity ratio. As explained previously, this ratio can take any value between 0 and 1. Since a β value of 1 defines a purely Newtonian fluid, we expect that fluids with viscosity ratios less than 1 should exhibit polymer behavior similar to increasing the Weissenberg number in polymer flow. The elastic vortex should therefore be larger for smaller β values because of its reliance on elasticity to form.
The elastic vortices in Figure 12 show the effect of decreasing β; from β = 0.4 to β = 0.1. The second elastic vortex, lower viscosity ratio, is larger than expected, indeed, significantly larger, roughly three times. Reducing β is predicted to have a similar effect to that of increasing the Weissenberg number, increasing the size of the elastic vortex. However, the presence of elasticity is not sufficient to produce elastic vortices; it was found that elastic vortices are only predicted to form for viscosity ratios not exceeding a certain value. In other words, a dilute polymer solution, mostly solvent with a viscosity ratio near 1, is always predicted to have Newtonian-like flow behavior without a vortex forming upstream of the bend corner. Conversely, a concentrated polymer solution, mostly polymer with a viscosity ratio near 0, is predicted to have secondary flows, vortices upstream of the bend corner, for a large range of Reynolds numbers. These predictions are in line with the experiments of Gulati [5] who found that stable elastic vortices form within the flow of semidilute DNA solutions in a 90 • micro bend channel above a certain threshold for the solution elasticity. Based on these findings we created a map to evaluate the effect of the three parameters, Wi, β, and El on the formation of elastic and inertial vortices. Figure 13 shows a phase diagram mapped for viscosity ratio versus elasticity number for Wi = 1.0. This value was chosen so the inverse of the elasticity number, El, is exactly the Reynolds number, Re. Based on these findings we created a map to evaluate the effect of the three parameters, Wi, β, and El on the formation of elastic and inertial vortices. Figure 13 shows a phase diagram mapped for viscosity ratio versus elasticity number for Wi = 1.0. This value was chosen so the inverse of the elasticity number, El, is exactly the Reynolds number, Re. Three regions form in this phase diagram: the 'devoid region', no vortex forms; the 'inertia region', downstream or inertial vortices form; and the 'elastic region', upstream or elastic vortices form. As can be seen in the diagram, the boundary between the elastic region and the devoid region is a horizontal line at about β = 0.6. This is interesting because it suggests that there is a critical viscosity ratio needed for the elastic vortex to form. It might be possible that this boundary is not asymptotic, and that it is dependent on the Weissenberg and the Reynolds numbers. Take note of the logarithmic x-axis. The effect of the Weissenberg number on this predicted flow map will be the subject of a follow-up paper.
chosen so the inverse of the elasticity number, El, is exactly the Reynolds number, Re. Three regions form in this phase diagram: the 'devoid region', no vortex forms; the 'inertia region', downstream or inertial vortices form; and the 'elastic region', upstream or elastic vortices form. As can be seen in the diagram, the boundary between the elastic region and the devoid region is a horizontal line at about β = 0.6. This is interesting because it suggests that there is a critical viscosity ratio needed for the elastic vortex to form. It might be possible that this boundary is not asymptotic, and that it is dependent on the Weissenberg and the Reynolds numbers. Take note of the logarithmic x-axis. The effect of the Weissenberg number on this predicted flow map will be the subject of a follow-up paper.
The graph in Figure 14 is the same phase diagram as Figure 13, but the x-axis is now the Reynolds number (the inverse of the elasticity number for Wi = 1.0). In the figure, it is easier to see that elastic vortices only form at small Reynolds number values and viscosity ratios less than 0.6. As it has already been shown that viscosity ratio has a similar effect on vortex development to that of increasing the Weissenberg number, the question arises as to whether or not β will delay the onset of inertial vortices. Looking at the figure, the boundary of the inertial region is essentially a vertical line (Re = 29-37) indicating that β is predicted to have only a very subtle effect on the onset of inertial vortices. In other words, the Reynolds number is predicted to be the parameter with a primary role in determining the onset of inertial vortices.  The graph in Figure 14 is the same phase diagram as Figure 13, but the x-axis is now the Reynolds number (the inverse of the elasticity number for Wi = 1.0). In the figure, it is easier to see that elastic vortices only form at small Reynolds number values and viscosity ratios less than 0.6. As it has already been shown that viscosity ratio has a similar effect on vortex development to that of increasing the Weissenberg number, the question arises as to whether or not β will delay the onset of inertial vortices. Looking at the figure, the boundary of the inertial region is essentially a vertical line (Re = 29-37) indicating that β is predicted to have only a very subtle effect on the onset of inertial vortices. In other words, the Reynolds number is predicted to be the parameter with a primary role in determining the onset of inertial vortices. Lastly, Figure 15 plots the vortex length of elastic vortices versus the elasticity numbers at which they form for three values of Weissenberg number (Wi = 0.5, 1.0, 1.5) and a viscosity ratio β = 0.1. The figure shows that for elasticity numbers between 0.1 and 10 the size of the elastic vortex is predicted to depend on both the elasticity number, El, and the Weissenberg number, Wi. On the other hand, the figure also shows that for El = 10 and above, the elastic vortex size is predicted to be of the same size for a given Weissenberg number independently of the elasticity number. These predictions indicate that the Weissenberg number plays the primary role in determining the size of the elastic vortex that forms upstream of the bend corner for a fixed viscosity ratio. As discussed earlier, as the Weissenberg number increases so does the size of the elastic vortex. These predictions are also consistent with the experiments of Gulati [5] who found that the Weissenberg number is the parameter that determines the presence and size of the elastic vortices forming within the flow of semidilute DNA solutions in a 90 • micro bend channel. Lastly, Figure 15 plots the vortex length of elastic vortices versus the elasticity numbers at which they form for three values of Weissenberg number (Wi = 0.5, 1.0, 1.5) and a viscosity ratio β = 0.1. The figure shows that for elasticity numbers between 0.1 and 10 the size of the elastic vortex is predicted to depend on both the elasticity number, El, and the Weissenberg number, Wi. On the other hand, the figure also shows that for El = 10 and above, the elastic vortex size is predicted to be of the same size for a given Weissenberg number independently of the elasticity number. These predictions indicate that the Weissenberg number plays the primary role in determining the size of the elastic vortex that forms upstream of the bend corner for a fixed viscosity ratio. As discussed earlier, as the Weissenberg number increases so does the size of the elastic vortex. These predictions are also consistent with the experiments of Gulati [5] who found that the Weissenberg number is the parameter that determines the presence and size of the elastic vortices forming within the flow of semidilute DNA solutions in a 90° micro bend channel.

Discussion
It has been shown during this research paper that higher polymer concentrations and higher inertia are both predicted to lead to larger vortices in a sharp bend geometry. Higher polymer concentrations lead to elastic vortices located upstream of the bend. Conversely, higher Reynolds numbers lead to inertial vortices located downstream of the bend.

Discussion
It has been shown during this research paper that higher polymer concentrations and higher inertia are both predicted to lead to larger vortices in a sharp bend geometry. Higher polymer concentrations lead to elastic vortices located upstream of the bend. Conversely, higher Reynolds numbers lead to inertial vortices located downstream of the bend.
It was also found that after vortices are formed their size is predicted to be primarily determined by the properties of the polymer solution and the Reynolds number. In the case of elastic vortices, the vortex size is set by both the Weissenberg number and the viscosity ratio. In the case of inertial vortices, the vortex size is primarily set by the Reynolds number with the Weissenberg number playing a secondary role.
What is more interesting is the fact that predictions indicate vortices can be suppressed by adjusting the properties of the polymer solution, elasticity and viscosity ratio, for a given solution volumetric flow rate. The transition to inertia driven vortices is predicted to occur somewhat abruptly at elasticity numbers equivalent to the critical Reynolds number for Newtonian flow. Above this critical Reynolds number, inertia vortices are predicted to form downstream of the bend corner with polymer solution properties affecting the size only moderately. Below this critical Reynolds number, all polymer solution properties are predicted to affect not only the size but the formation of elastic vortices upstream of the bend corner. The size of the vortex is predicted to be determined by the Weissenberg number while the formation of the vortex itself is predicted to be determined by the viscosity ratio, i.e., the polymer concentration in the solution.
It should be noted that with our numerical approach we did not expect to predict or capture any secondary flow that might be induced or be present in the actual 3D geometry used as reference [5]. Similar works have found or predicted a minor effect of the aspect ratio of the channel in the size of the lip vortex [21] or the presence of secondary flows [28].
As the present work summarizes, our current primary interest was to evaluate the role that elasticity plays in the formation of lip vortices in the flow around a sharp corner. It is clear that modifying the underlying properties of the polymer solution predicts that the fluid dynamics can be significantly altered. Vortices, elastic and inertial, could both be suppressed or promoted by adjusting the polymer solution properties for a given volumetric flow situation.

Further Work
The phase diagrams presented in Figures 13 and 14 were created for a single Weissenberg number. The non-dimensional elastic vortex length behavior presented in Figure 15 was created for one single viscosity ratio. Additional diagrams for different values of both might confirm the universality of the predicted behavior.
In the present work, we used the FENE-P model to predict the behavior of the polymer solution. This model is computationally stiff. Extreme care has to be taken to find fully converged numerical solutions. Replicating the work with other polymer models would further verify the present findings.
In the long term, our goal is exploring the possibility of developing polymer and polymerlike solutions that respond to and change properties in a controlled manner based on self-induced flow instabilities. Ideally, developing an experimental setup where polymer flow could be tested to create a real phase diagram for the fluid dynamics would validate what the current work promises, controlling the formation of vortices based on the polymer solution's properties.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
As discussed within the text, a 10:1 branches-to-channel height ratio for the geometry under analysis was found to be sufficient to predict fully developed flow before the bend corner. Figure A1 presents a comparison of the streamwise velocity component at two different locations. We are comparing profiles at a distance equivalent to 7H and 8H from the inlet. Figure A1a presents the comparison for the same elastic vortex case discussed in the grid convergence analysis outlined in Figure 3; recall that an elastic vortex forms before the bend corner. The two profiles are nearly identical, indeed indistinguishable. The weighted average difference between the two profiles is 0.02%. Figure A1b presents the comparison for the same inertial vortex case discussed in the grid convergence analysis section; recall that an inertial vortex forms after the bend corner. The two profiles are also nearly identical. The weighted average difference between the two profiles is 0.52%. grid convergence analysis outlined in Figure 3; recall that an elastic vortex forms before the bend corner. The two profiles are nearly identical, indeed indistinguishable. The weighted average difference between the two profiles is 0.02%. Figure A1b presents the comparison for the same inertial vortex case discussed in the grid convergence analysis section; recall that an inertial vortex forms after the bend corner. The two profiles are also nearly identical. The weighted average difference between the two profiles is 0.52%. In both cases, it is clear the chosen length to height ratio is sufficient to obtain the fully developed flow condition before the bend corner. Hence, the inlet boundary condition should have no influence or effect on the size predictions of both elastic and inertial vortices.
As in figure A1, figure A2 presents a comparison of the streamwise velocity component at two different locations. We are now comparing profiles at a distance equivalent to 7H and 8H after the bend corner, or 3H and 2H far from the geometry outlet, respectively. Figure A2a presents the comparison for the same elastic vortex case discussed above. The two profiles are nearly identical, again indistinguishable. The weighted average difference between the two profiles is 0.01%. Figure A2b presents the comparison for the same inertial vortex case discussed above. The two profiles are again nearly identical. The weighted average difference between the two profiles is 0.65%. In both cases, it is clear the chosen length to height ratio is sufficient to obtain the fully developed flow condition before the bend corner. Hence, the inlet boundary condition should have no influence or effect on the size predictions of both elastic and inertial vortices.
As in Figure A1, Figure A2 presents a comparison of the streamwise velocity component at two different locations. We are now comparing profiles at a distance equivalent to 7H and 8H after the bend corner, or 3H and 2H far from the geometry outlet, respectively. Figure A2a presents the comparison for the same elastic vortex case discussed above. The two profiles are nearly identical, again indistinguishable. The weighted average difference between the two profiles is 0.01%. Figure A2b presents the comparison for the same inertial vortex case discussed above. The two profiles are again nearly identical. The weighted average difference between the two profiles is 0.65%. In both cases, it is clear the chosen length to height ratio is enough to predict that the flow will recover completely from the perturbation caused by the vortices before the geometry outlet, nearly reaching the fully developed flow condition. Hence, the outlet boundary condition should also have no influence or effect on the size predictions of either elastic or inertial vortices. In both cases, it is clear the chosen length to height ratio is enough to predict that the flow will recover completely from the perturbation caused by the vortices before the geometry outlet, nearly reaching the fully developed flow condition. Hence, the outlet boundary condition should also have no influence or effect on the size predictions of either elastic or inertial vortices.

Appendix B
As discussed within the text, the size (length) of all predicted vortices was determined numerically by the change in direction of the tangential velocity near the wall. We indicated that doing this was equivalent to a change in sign of the shear stress or skin friction. Figure A3 presents a comparison of the wall shear stress and the near wall streamwise (tangential) velocity. To allow a direct comparison, we normalized both stress and velocity by the maximum value along the respective walls. Figure A3a presents the comparison for the same elastic vortex case discussed in Appendix A. The two profiles cross the zero line at nearly the same x-coordinate location. The length of the vortex determined by either method is the same with a 0.0% deviation. Figure A3b presents the comparison for the same inertial vortex case discussed in Appendix A. The two profiles cross the zero line at almost the same y-coordinate location. The length of the vortex determined by either method is roughly the same, with a 0.5% deviation. In both cases, elastic and inertial vortices, it is clear that determining the vortex size by the change in direction of the tangential velocity near the wall is equivalent to a change in the sign of the wall shear stress. Any deviation between the methods that could be predicted is extremely small compared to the actual vortex sizes that were calculated. In both cases, elastic and inertial vortices, it is clear that determining the vortex size by the change in direction of the tangential velocity near the wall is equivalent to a change in the sign of the wall shear stress. Any deviation between the methods that could be predicted is extremely small compared to the actual vortex sizes that were calculated.