The Viscoelastic Swirled Flow in the Confusor

A two-dimensional mathematical model for a steady viscoelastic laminar flow in a confusor was developed under the condition of swirled flow imposed at the inlet. Low density polyethylene was considered as a working fluid. Its behavior was described by a two-mode Giesekus model. The proposed mathematical model was tested by comparing it with some special cases presented in the literature. Additionally, we propose a system of equations to find the nonlinear parameters of the multimode Giesekus model (mobility factor) based on experimental measurement. The obtained numerical results showed that in a confusor with the contraction rate of 4:1, an increase in the swirl intensity at Wi < 5.1 affects only the circumferential velocity, while the axial and radial velocities remain constant. The distribution pattern of the first normal stress difference in the confusor is qualitatively similar to the one in a channel with abrupt contraction, i.e., as the viscoelastic fluid flows in the confusor, the value of N1 increases and reaches a maximum at the end of the confusor. Dimensionless damping coefficients of swirl are used to estimate the swirl intensity. The results show that the swirl intensity decreases exponentially.


Introduction
Investigation of viscoelastic fluid flows in channels of various configurations is of particular practical interest for polymer production. During extrusion, the polymer material passes through a channel with a screw, thus creating a swirling flow directed to the die nozzle. The nozzle geometry varies depending on the shape of the finished product, but a convergent channel (confusor) is an integral part of die [1]. Since we only consider the flow of initially swirled viscoelastic fluid in the convergent channel, our review is limited to the results covering this specific subject.
The literature mainly considers individual extreme cases, namely viscoelastic flows in channels with abrupt contraction (planar and axial symmetry flow), flow patterns in viscoelastic fluids flowing in a limited space with a rotating wall (bottom or upper lid, pipe surface), and cone-and-plate flows [2].
Interestingly, the research of viscoelastic fluid flows in channels with abrupt contraction is mainly concerned with planar flows, while the planar configuration is better suited to visualization studies through birefringence strand techniques and particle image velocimetry (PIV) [3,4]. Flow in a channel with abrupt contraction (as well as flow past a cylinder [5]) is a well-known benchmark problem used for the testing of reliability of new or modified numerical methods simulating the flows of viscoelastic fluid that employ multimodal rheological equations of state, e.g., the Giesekus model [6], the Phan-Thien-Tanner model (PTT) [7], Extended Pom-Pom [8], Oldroyd-B [9], and Fene-P [10]. A special feature of viscoelastic flows in channels with abrupt contraction is a recirculation zone, the size of which depends on the Weissenberg/Deborah number [11] and on extensional-flow properties [12]. The comprehensive review presented in [13] summarizes key factors influencing secondary entry flows for polymer melts, while outstanding issues in numerical methods and novel and challenging applications of viscoelastic fluids are discussed in [14].
Earlier [19], it was found that the planar contraction flow gives rise to very low vortex activity in the salient corner, unlike similar flows in circular contractions. Therefore, a circular channel is the most interesting geometry, particularly as far as swirled flows are concerned, because both the geometry and the system of equations describing hydromechanical processes are invariant to a circumferential angle phi (angle of rotation about a symmetry axis in a cylindrical coordinate system).
According to the literature, viscoelastic fluid flows in channels with abrupt contraction have been studied extensively for the case when a developed velocity profile (Newtonian or non-Newtonian) is set at the channel inlet, and when the inlet and outlet parts of the channel are long enough so as not to affect the flow pattern near the contraction.
Lately, embedded software (e.g., ANSYS-Polymat) has become an increasingly popular tool for the estimation of both the discrete spectrum of relaxation and non-linear parameters of differential PTT and Giesekus models [1,29]. According to preliminary analysis, our proprietary software approximates the viscosity curve with fewer modes and the same error of approximation of numerical and experimental data. Normal stresses exist in viscoelastic media, as discovered in experiments by Garner and Nissan [30] and interpreted by Weissenberg [31]. These results triggered the research of the structure of swirled viscoelastic flows in a confined cylinder with a rotating bottom lid [32][33][34][35][36][37]) or pipe wall [38]. Laminar pipe flow with a controllable wall swirl has been studied in [39] to explore the behavior of inelastic shear-dependent fluids. A vortex shedding regime was illustrated using experimental data in [35]. It was also observed that the dimensionless circumferential velocity decreases with the increase in the Weissenberg number, We [34]. The obtained results were employed for the development of advanced rotary rheometers with plate-plate and cone-plate measurement systems.
Three-dimensional numerical simulation [33], contrary to earlier experiments in [32], demonstrated that the structure of swirled flow in a confined cylinder with a rotating bottom lid is axisymmetric. It should be mentioned that numerical results obtained in [33] were validated by checking the stability criterion [40] for the case of highly unsteady spiral vortex flow of viscoelastic fluid. Thus, our two-dimensional approach to the construction of a mathematical model for a swirled flow of two-mode Giesekus fluid is consistent with the physical pattern of flow.
Numerical analysis [1] of screw swirling effects on fiber orientation in large area additive manufacturing polymer composite deposition is the most similar in its content to our study. The authors [1] considered a 3D problem using the exponential form of the Phan-Thien-Tanner model (PTT) and commercial software ANSYS PolyFlow. They had to use three-dimensional formulation because it is impossible to modify the standard rheological equation of state in ANSYS PolyFlow when using an axisymmetric problem statement to take fluid rotation into account.
The present work aims to develop a two-dimensional mathematical model of swirled viscoelastic fluid flow in a circular convergent channel (confusor). Using this model, the distribution of hydrodynamic parameters is obtained more easily, and it complies with the results obtained for a 3D problem.

Materials and Methods
Let us consider a steady-state swirled flow of a two-mode viscoelastic fluid in a confusor with a contraction rate of 4:1 (Figure 1). At the inlet, the swirled flow is prescribed by the boundary conditions, which are invariant with respect to variable ϕ, so the main problem can be reduced to a two-dimensional form since the governing equations for the considered problem are also invariant to variable ϕ. Thereby, the mathematical model of viscoelastic fluid flow in the confusor is as follows: ∂V r ∂r where V r , V ϕ , V z are radial, circumferential, and axial velocity, respectively; r, ϕ, z are cylindrical coordinate system variables (axis z is the axis of rotation); p is pressure, Boundary conditions: The shear stresses and pressure are assumed to be zero at the outlet.
Here, V a = Q/ πR 2 1 is a mean velocity over the channel cross-section, and Q is the flow rate (m 3 /s).
Boundary condition (5) is an ideal model with K = ωR 1 /V a . In this study, we simplified this condition as follows: In this study, we use the two-mode Giesekus model [7].
For the case of stationary flows ∂σ ∂t = 0, then the upper convective derivative takes the form is the strain rate tensor; (λ m ,η m ) are the relaxation spectra, α m is the rheological parameter of the Giesekus model.
Equation (8), written in a cylindrical coordinate system for the considered case (axisymmetric formulation), is presented in the Appendix A.
As a specific liquid, we are going to consider DSM Stamylan LD 2008 XC43 lowdensity polyethylene (LDPE) melt with ρ f = 921 [kg/m 3 ] [6]. The authors of [5] defined the parameters of the four mode Giesekus model, but such a high number of modes may cause a convergence problem. So, we defined the parameters of the two-mode Giesekus model (Table 1) for this fluid using the following algorithm. The pairs of λ i , η V i and η N were found from:  We used the relative deviations of the experimental and calculated values of the dynamic moduli that allowed the improvement of the accuracy of the approximation [41]. The parameters α k (k = 1, m) were calculated from: where . . , n. j = 1...n is the number of an experimental point, and k = 1...m is the mode number. The relation between the shear stresses of the k-th mode and the shear rate for a torsional flow of Giesekus fluid between two parallel plates (measurement system of rheometer) can be written as follows: (∀k = 1, . . . , m and j = 1, . . . , n), γ is the shear rate (1/s), and τ is the shear stress (Pa).

Approbation
The obtained parameters of the two-mode Giesekus model were tested on a classic problem of a viscoelastic fluid flowing around a cylinder located between two infinite plates [5]. Numerical results were obtained using the OpenFoam with the viscoelastic package (planar flow). As Figure 3 shows, the distribution and value of the velocity components in the flow according to the two-mode Giesekus model are consistent with both experimental data and numerical results obtained using the four-mode Giesekus model. The numerical implementation of the problem (1)-(8) is carried out in the Comsol Multiphysics package, which allows us to solve the custom equations by using the partial differential equations (PDE) package. The computational domain of the channel was subdivided using quadrangular elements with a minimum element quality of 0.78; the total number of elements was 45,300. The PARDISO method was used as a solver. The problem was solved on the XeonGold computational server with 24 cores and 512 Gb RAM.
As an approbation of our mathematical model (1)

Results and Discussion
In this work, we considered the fixed geometry of the confusor: Figure 1). According to the preliminary data, the length of the outlet part is L 3 = 10D 2 and it is sufficient for the outlet boundary conditions to have no effect on the solution inside the confusor. Figure 5 shows the profiles of the normalized velocities: axial (a), circumferential (b), and radial velocity (c). The profiles shown in the figure are plotted in cross-sections of the channel (Figure 1), equidistant from each other at a distance of 0.01 m. It should be mentioned that the z-axis of the cylindrical coordinate system is centered in the crosssection corresponding to the outlet of the confusor and oriented in the flow direction. The components of the velocity vector are normalized by the mean velocity calculated for each cross-section V a (i) (i indicates the cross-section of the confusor). It is known that at a constant liquid flow rate, a decrease in the channel cross-section leads to an increase in mean velocity. For example, if for z = −0.045 (inlet) the velocity is given as V a (inlet) = 0.03 (m/s), then at z = −0.02 (R 1 ≈ 0.01248 (m)) the mean velocity will increase to V a (z=−0.02) ≈ 0.07704 (m/s). All primary data were processed in Excel to reduce the accumulation of errors. Figure 5 shows a comparison of the obtained results (K = 6) with similar results in a confusor without swirling (K = 0). As can be seen from the figure, the profile of the normalized axial velocity (Figure 5a) stretches as the fluid flows in the confusor due to the channel narrowing, then drops sharply and tends to the distribution corresponding to the steady flow in a round pipe. The behavior of axial velocity is similar to that in a channel with sudden contraction. An analysis of the obtained data revealed that the swirl intensity (parameter K in Equation (7)) affects only the circumferential velocity (V ϕ ), the distribution of which sharply tends to zero as the fluid flows in the confusor, i.e., the swirl intensity drops abruptly. The presence of radial velocity is a characteristic feature of fluid flow in a convergent channel (confusor) (Figure 5c). The figure demonstrates that the radial velocity is comparable to the axial one up to the value of z = −0.01 (3/4 of the confusor length). Normal stress profiles for viscoelastic fluid flow in a convergent channel differ significantly from the similar profiles for fluid flow in a straight pipe (Figure 5c,d). For example, as the fluid flows in the confusor, the component T zz (a component of the extra stress tensor) takes on progressively higher values, with the maximum located in the middle region between the axis of symmetry and the channel wall. The T rr component over the channel cross-section, on the contrary, takes negative values, the highest absolute value of which is also concentrated in the region between the axis of symmetry and the channel wall. The presence of flow swirling (K = 6) insignificantly increases the value of the local extremum, which is consistent with the distribution of the velocity components.
The following figure illustrates the influence of the Weissenberg number on the distribution of axial velocity and normal stresses when the swirling viscoelastic fluid flows in a confusor. We considered Wi = 1.68 and Wi = 5.04 that were calculated by the formula Wi = λV a (R 2 ) −1 for the velocities V a (inlet) = 0.01 m/s and V a (inlet) = 0.03 m/s, respectively. Here, the value V a (inlet) is calculated for R 1 (inlet) = 0.02 (m). In the considered cross-sections, the greatest difference in the profiles of the normalized axial velocity is observed at z = −0.03, i.e., in the inlet region, which is apparently related to the selected boundary conditions at the channel inlet (Figure 6a). It was found (Figure 6b,c) that with an increase in the Weissenberg numbers, there is a significant increase in normal stresses, especially in the outlet region of the confusor. In this case, the stress profile with an increase in the Weissenberg number is characterized by the presence of a pronounced local extremum: for T zz this is a local maximum, for T rr this is a local minimum. In particular, for z = −0.01 the ratio T zz (max) /T zz (axis symmetry) = 1.47 for Wi = 5.04 and |T rr | (max) /|T rr | (axis symmetry) = 1.12 for Wi = 1.68.  It should be noted that, here, the current values of the axial velocity were also normalized by the actual mean velocity corresponding to a considered cross-section (meaning that as the fluid flows in the confusor, the mean integral flow velocity for the entire section increases due to the channel narrowing). As can be seen from Figure 7a, the flow of viscoelastic fluid is characterized by an increase in the axial velocity on the channel axis as it flows in the confusor, then a sharp drop and stabilization of the flow, i.e., after the outflow from the confusor, the velocity distribution tends to the one in steady pipe flow. The obtained results are in good agreement with similar results in the channel with sudden contraction. As can be seen from the figure, an increase in the Weissenberg number leads to a sharper increase in the axial velocity inside the confusor and the presence of a local maximum in the output region. As shown by the numerical results, for the considered case, the presence of flow swirl does not affect the distributions of both the axial velocity and the first normal stress difference. Note that N 1 (i−max) was calculated for each distribution individually, so the maximum value of N 1 /N 1 (i−max) does not exceed unity. Similar to the axial velocity distribution, an increase in the Weissenberg number leads to a sharper increase in the value of N 1 /N 1 (i−max) . When studying swirling flows, it is convenient to use Equation (12) to estimate the swirl intensity [43]. where It should be noted that when substituting v ϕ = V ϕ /V a and v z = V z /V a in Equation (9), the value V a cancels out, so there is no need to recalculate V a for the cross-sections inside the confusor. For the convenience of analysis, Figure 8 shows the distribution of m(z)/K, which made it possible to bring together the investigated relationships for K = 2, 4, 6. It can be seen from the figure that the dependences m * (z) = m(z)/K for Wi = 5.04 and Wi = 1.68 are satisfactorily described by the formula of the form m * (z) = Ae B , which is consistent with the results of [43], in which the damping of the power-law fluid flow in a circular tube was studied. The calculations have shown that the presence of a confusor leads to a faster decrease in the swirl intensity. For example, in Figure 8, the results for the similar fluid flow in a straight pipe with the same boundary conditions at the inlet (Wi = 5.04) are shown in green. Thus, the presence of a convergent channel leads to the suppression of the swirl intensity and the flattening of the velocity profile.

Conclusions
In this paper, we developed a two-dimensional mathematical model of a steady laminar viscoelastic flow in the confusor under the condition of swirled flow imposed at the inlet. The parameters of the two-mode Giesekus model were obtained using the proposed relation. Approbation of the mathematical model showed good agreement with the literature data. It was observed that swirl intensity imposed at the inlet of the confusor with a contraction rate of 4:1 does not affect axial and radial velocities. The distribution of axial velocity along the axial direction is the following: it stretches along the channel axis and then becomes more flattened in the outlet section of the confusor. The magnitude of radial velocity is comparable to the magnitude of axial velocity up to 3/4 of the confusor length. The circumferential velocity distribution depends on boundary conditions; it decreases extensively along the axial direction and nearly disappears after 1/2 of the confusor length. A significant increase in normal stresses is observed with an increase in the Weissenberg number. Near the outlet of the confusor, the magnitudes of the T zz component in the central region prevail over the corresponding magnitudes in the near-boundary region, which is fundamentally different from the flow of viscoelastic fluid in a straight pipe. The post processing of numerical results showed that the damping of the swirled flow in the confusor is more intense compared to a straight pipe. In this case, the distribution of damping coefficients of swirling along the length of the confusor can be described by an exponential function. The proposed mathematical model makes the estimation of hydrodynamic parameters a lot easier compared with a 3D statement.
Author Contributions: Methodology, data analysis and curation, A.K.; numerical simulation R.Z.; literature review, J.K.; supervision E.V. All authors have read and agreed to the published version of the manuscript.
Funding: This study was supported by the Russian Science Foundation (Project no. 19-11-00220).
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.

Acknowledgments:
The authors finally want to thank the reviewers for informative comments and advice that did improve the paper markedly.

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