An Effective Interface Tracking Method for Simulating the Extrudate Swell Phenomenon

The extrudate swell, i.e., the geometrical modifications that take place when the flowing material leaves the confined flow inside a channel and moves freely without the restrictions that are promoted by the walls, is a relevant phenomenon in several polymer processing techniques. For instance, in profile extrusion, the extrudate cross-section is subjected to a number of distortions that are motivated by the swell, which are very difficult to anticipate, especially for complex geometries. As happens in many industrial processes, numerical modelling might provide useful information to support design tasks, i.e., to allow for identifying the best strategy to compensate the changes promoted by the extrudate swell. This study reports the development of an improved interface tracking algorithm that employs the least-squares volume-to-point interpolation method for the grid movement. The formulation is enriched further with the consistent second-order time-accurate non-iterative Pressure-Implicit with Splitting of Operators (PISO) algorithm, which allows for efficiently simulating free-surface flows. The accuracy and robustness of the proposed solver is illustrated through the simulation of the steady planar and asymmetric extrudate swell flows of Newtonian fluids. The role of inertia on the extrudate swell is studied, and the results that are obtained with the newly improved solver show good agreement with reference data that are found in the scientific literature.


Introduction
Free-surface flows are encountered in many polymer processing and environmental applications [1,2]. Nevertheless, the variety of the analytical solutions for free-surface flows is usually very limited, even for very simple cases [3]. On the other hand, the experimental observations of real phenomena are onerous [3], and many experimental techniques are suitable for single-phase flows and undergo many difficulties to be extended to two-phase flows [4]. For these reasons, the use of numerical simulations to provide useful information about free-surface flows would be of great advantage. Flows with a free-surface are difficult to be modeled since the free-surface is a moving boundary, whose location is merely known initially, and it has to be determined later during the simulation [5]. There are different ways of modeling free-surface flows: the Interface Tracking (IT) approach, in which the free-surface is tracked using a sharp interface, and a dynamic computational grid is applied to follow the movement of the free-surface; and, the Interface Capturing (IC) approach, in which the free-surface is not treated as a sharp interface and, generally, the computational grid is static. Among the methods following the IC approach, the Marker-And-Cell (MAC) [6] method is based on a finite difference scheme applied to an Eulerian grid to solve the Navier-Stokes equations for the fluid flow die temperature is lower than the melted polymer, the viscosity of the melted polymer increases at the wall and, therefore, the flow becomes limited and it undergoes a lower swell when compared to the case where the die temperature is higher than the melted polymer, where the liquid is less viscous at the wall, which lubricates the flow and results in a larger swell. Another important parameter for the polymer extrudate swell is the molecular structure of the polymer, since the first normal stress difference that appears at the die exit is affected by the molecular weight distribution. Furthermore, long chain branching enhances extrudate swell in polymers [33].
It is crucial to know the dimension of the emerging polymer from the die, since it determines the exact size of the extruded products; therefore, many attempts have been made to obtain equations for anticipating the swell ratio [31,34]. As an example, Tanner [35] introduced an equation for the swell ratio of a Maxwell-type constitutive equation of a viscoelastic fluid, and then revised it several years later [36] while considering that the first normal stress difference obeys N 1 = kτ m , where τ is the fluid shear stress and k and m are model parameters, instead of the initial assumption N 1 = constant × τ 2 . A large number of investigations of the die swell phenomenon are reported in the literature, which goes from two-dimensional (2D) to three-dimensional (3D) flows, from Newtonian to viscoelastic fluids rheology or from isothermal to non-isothermal flows. Crochet and Keunings [21] presented a 2D numerical investigation on slit die swell with three finite element meshes and two different techniques, a mixed and an extended u − v − p methods, and concluded that the swelling ratio depends upon the method that is used for its calculation and, for a given method, it is highly dependent upon mesh refinement. Subsequently, Mitsoulis et al. [22] used viscometric flow equations to simulate the swelling of viscoelastic fluids from long slit and capillary dies. The obtained results showed that, even when using this simple model for viscoelastic calculations, they were in good qualitative agreement with other numerical simulations, which employ a constitutive equation satisfying tensorial invariance, but accelerate the breakdown of the iterative scheme. Another work from the group of J. Vlachopoulos was presented by Karagiannis et al. [23], which studied the 3D free surface die swell of a Newtonian fluid in different geometries, specifically square, rectangular, equilateral triangular, bow-tie, and key-hole-shaped geometries. The obtained results were compared with experimental measurements and other numerical calculations with favourable agreement. The swelling ratios were found to strongly depend on the die geometry. The group of Crouchet also presented calculations for steady state 3D free surfaces of Newtonian and power-law fluids in the work of Wambersie and Crochet [24]. They combined a pseudo-transient marching technique, a decoupling algorithm, and a conjugate gradient solver to reduce the cost of the 3D calculations. The method was employed to study the circular, square, and rectangular die-swell problems, where the effects of inertia and shear-thinning were revealed. The works of Legat and Marchal [25,26] addressed the prediction of 3D free surface extrudate flows with a fully implicit finite element algorithm, in the sense that a Newton-Raphson scheme was applied to all variables and is geometrically general. The algorithm was employed to compute the extrudate swell of a rectangular die and in various complex sections containing multiple corners. The obtained results showed that the extrudate shape exhibits large deformations in the vicinity of all re-entrant corners, which would not be possible to predict in 2D simulations. Subsequently, the works of Georgiou and Boudouvis [27] and Mitsoulis et al. [28] studied the effects of inertia, surface tension, gravity, slip, and compressibility for both the 2D planar and axisymmetric extrudate-swell flows of Newtonian fluids. Recently, 3D isothermal and non-isothermal viscoelastic flow calculations with a transient finite element method for predicting extrudate swell of domains containing sharp edges were conducted by Spanjaards et al. [29,30]. The obtained results showed that the extrudate swelling is highly dependent on the rheological parameters and the constitutive model used, and that the wall temperature of the die can lead to a change in the bending direction. All of the presented works employed the finite element method to solve the problem of the extrudate swell for Newtonian and viscoelastic fluids. The present work aims to revisit the Newtonian extrudate swell flow problem by using the finite volume method, which is the core of the open-source computational library OpenFOAM [37].
When dealing with a steady state process, which is the case of profile extrusion, IT is usually the best alternative, since it does not present the problems that are related to interface diffusion inherent to the IC methods. OpenFOAM [37] comprises a solver to simulate free-surface flows following an IT approach, which was proposed by Tuković and Jasak [38]. One of the disadvantages generally raised to the IT methods is their computational efficiency, and many attempts have been carried out in order to increase the convergence rate of the pressure-velocity calculations [39][40][41]. Tuković et al. [41] proposed a non-iterative Pressure-Implicit with Splitting of Operators (PISO) algorithm based on extrapolation of mass flux, nodal velocity, and pressure from two previous time steps in order to have an approximation of these quantities in the new time step, and obtained a second order temporal accuracy in the cases with static and dynamic meshes.
This work aims to assess the capability of the solver that was developed by Tuković and Jasak [38] with the non-iterative PISO algorithm proposed by Tuković et al. [41] to efficiently simulate the extrudate swell phenomenon. For this purpose, the developed solver couples the interfaceTrackingFvMesh and interTrackMeshMotion libraries that are available in OpenFOAM [37] with the consistent second-order time-accurate noniterative PISO algorithm. A least-squares volume-to-point interpolation method for the grid movement, which enables an efficient and accurate tracking of the free-surface motion, is employed. The enhanced algorithm is used to simulate the steady-state Newtonian extrudate swell problem in both planar and axisymmetric geometries for a parametric study of the effects of inertia, and the obtained results are compared with the reference data of Mitsoulis et al. [28]. Notice that, although the results presented here are limited to Newtonian fluids, the effects which are discussed can be qualitatively applied to all fluids (e.g., viscoelastic fluids). The main aim of this work is to present an open source finite volume numerical framework that can handle the extrudate swell problem in an efficient way, which can be extended in the future to allow simulating other fluids rheology.

Mathematical Formulation
In this section, the governing equations and numerical method employed for simulating the two-phase fluid flow with a sharp interface are described, while using the finite-volume method and an interface tracking algorithm for the moving mesh. The numerical scheme developed in this work is enhanced with the consistent second-order time-accurate non-iterative PISO algorithm [39][40][41] to reduce the computational wall time of the simulations and, for the moving mesh calculation, a Laplacian scheme is used with a least-squares interpolation, which allows for robust and stable deformation of the interface.

Governing Equations
The mass and linear momentum conservation laws are the equations of motion governing the isothermal flow of incompressible Newtonian fluids inside an arbitrary volume V bounded by a closed surface S, where n is the outward pointing unit normal vector on S, u is the fluid velocity, u s is the surface S velocity, ν is the fluid kinematic viscosity, and P is the kinematic pressure obtained by subtracting the hydrostatic kinematic pressure P hydrostatic = g · r from the absolute pressure, where g is the gravitational acceleration and r is the position vector.
For an arbitrary moving volume, the relationship between the rate of change of the volume V and the velocity u s is defined by the geometrical (space) conservation law [42], When considering that the fluid phases are immiscible, the fluid flow Equations (1) and (2) can be used for each phase individually and, at the interface, the proper boundary conditions must be used.

Kinematic Condition
The kinematic condition states that the fluid velocities on the two sides of the interface, u 1 f and u 2 f , must be continuous (see Figure 1),

Dynamic Condition
From the momentum conservation law follows the dynamic condition, which states that forces acting on the fluid at the interface are in equilibrium. The general form of the dynamic condition at the interface, which gives the fundamental relationship between the jump in stress across an interface and the surface tension force, is given by, where T 1 and T 2 are the stress tensors that are defined in terms of the local fluid pressure and velocity fields, as respectively, σ is the interfacial tension and ∇ s = [I − nn] · ∇ = ∇ − n ∂ ∂n is the tangential gradient operator, which appears because σ and n are only defined on the surface. Equation (5) is a vectorial equation, which is often written in terms of its normal and tangential components. We proceed by deriving the normal and tangential force balances appropriately at a fluid-fluid interface that is characterized by an interfacial tension σ.
From the normal force balance follows the pressure jump across the interface [43], where κ = −∇ s · n is twice the mean curvature of the interface. The tangential force balance yields a relation between the normal derivative of tangential velocity on the two sides of the interface [43], where u t = (I − nn) · u is the tangential velocity component.

Numerical Method
The numerical integration in time of the mathematical model that is described in Section 2.1 is performed using a second order accurate implicit method, and the integral form of the fluid flow equations are discretized in space using a second order accurate cell-centred unstructured finite volume method. Detailed information of the computational domain discretization, of the mathematical model, and interface tracking procedure can be found in Tuković and Jasak [38], being out of the scope of this work. Here, we devote our attention to the improvements performed in the numerical algorithm related to both efficiency and robustness of the calculations. For this purpose, the consistent non-iterative PISO algorithm [39][40][41] was employed to assure the pressure-velocity coupling in the calculation of the free surface flow studied in this work, and a least-squares volume-topoint interpolation method to compute the interface motion, were newly-implemented in the two-phase fluid-flow solver with a sharp interface.

Consistent Non-Iterative PISO Algorithm
The rate of convergence of the collocated PISO algorithm is known to be problem dependent. In this segregated algorithm, a velocity correction term is neglected, which affects the path to convergence or may either cause the divergence of the numerical simulation, as a result of an exaggerated pressure correction field. Nevertheless, a common remedy for alleviating this problem is to under relax the pressure field. However, the rate of convergence remains a problem. In this work, we modified the original twophase flow interface tracking solver, which is based in the segregated PISO algorithm, by approximating the velocity correction at the main grid point by a weighted average of the velocity corrections at the neighboring locations, the so-called consistent PISO algorithm [39][40][41].
A brief summary of the collocated PISO algorithm is described hereafter, along with the modification performed to the original formulation, by using the consistent counterpart of the algorithm, to improve the convergence rate of the calculations. A detailed analysis of both algorithms can be found in Van Doormaal et al. [39], Issa [40] and Tuković et al. [41]. First, the discretized momentum equations are given by with which is a weighted average that consists of the contribution of the neighbor cell with centroid F, a u F , and the contribution of the current cell with centroid C, a u C , to the velocity of the neighbour cells u F . Notice that f refers to a face of the current cell, which shares it with a neighbor cell. The transient term contribution D u C and the source term contribution B u C are defined as the vector operators: where V C is the volume of current cell C. Equation (8) is solved to obtain a momentum conserving velocity field u * . Subsequently, the mass flow rate,ṁ * f , at the computational element faces should be updated using the Rhie-Chow [44] interpolation, which allows for obtaining a momentum satisfying the mass flow rateṁ * . Here, S f denotes the normal face area vector and all of the values with an over bar are obtained by linear interpolation between the values at points C and F. Subsequently, we assemble the pressure correction equation [40], and solve it to obtain a pressure correction field P . In the original PISO algorithm, the last term in the RHS of Equation (13) is neglected, which affects the convergence rate, because, the larger this term, the higher the error present in the approximation at each iteration is. Finally, the mass flow rate at the element faces (ṁ * * f ) and the pressure (P * C ) and velocity (u * * C ) at the element centroids are updated with the corrected pressure field P by, where the superscript (n) denotes the solution at time t = n and α p is the under relaxation factor for the pressure correction values, which increases the robustness and convergence behavior of the PISO algorithm. Nevertheless, even when using under-relaxation factors in the PISO algorithm, the convergence rate is problem dependent. To improve the efficiency of the two-phase fluid flow calculations, we modified the original PISO algorithm by simply assuming that the velocity correction at point C is the weighted average of the corrections at the neighboring points, which can be written as, Hence, the neglected term in the PISO algorithm (last term in the RHS of Equation (13)) is replaced by the approximate value that is given by Equation (18). Thus, in the consistent PISO algorithm a smaller term is discarded, which allows for obtaining more accurate velocity corrections with the momentum equations. Therefore, the convergence rate of the consistent PISO algorithm is higher than the one of the original PISO algorithm. Notice that our approach to simulate steady state free surfaces is, in fact, a time-dependent marching technique that allows separately calculating the free surface movement and the other fields at the different time steps, which reduces the number of iterations needed to converge to a steady solution [24].
A detailed analysis regarding the numerical setup for the PISO algorithm allowed for concluding that 10 outer corrector loops and three pressure-velocity correctors were needed to obtain stable and converged iterative solutions for all the cases simulated in this work. Additionally, for the PISO algorithm, the simulations only converged at maximum with a Courant number of 0.2, while using the consistent PISO algorithm the simulations were performed with a Courant number of 1.

Least-Squares Volume-to-Point Interpolation
The mesh deformation is calculated using the Laplace mesh motion equation with variable diffusivity [45]. The method discretizes the motion equation using the cell-centred FV method, by which vertex displacements are obtained using a reconstruction procedure, instead of the commonly employed FE discretization. Detailed information regarding this procedure can be found in Tuković and Jasak [38] and Jasak and Tuković [45].
Following the work of Tuković et al. [46], we employ the weighted least-squares method and linear fitting function to reconstruct the vertex displacements from the cellcentre displacements of the cells surrounding the vertex, which increased the robustness of the free-surface flows calculation. It is worth noting that, without the least-squares volumeto-point interpolation, the original mesh motion algorithm available in OpenFOAM [37] was not able to compute the mesh deformation that occurs in the simulations of the extrudate swell flows presented in Section 3.
Consider the interpolation stencils for the internal (i) and the boundary (b) vertices that are given in Figure 2. The former is constituted of all cells sharing the vertex, while, in the latter, the boundary faces sharing the corresponding vertex are also included into the stencil. In the vicinity of each vertex i (or boundary b), a linear interpolation function is considered: where C i is the unknown coefficient vector, and the field value φ i0 and the reference position r i0 are obtained as the weighted average of cell-centre field and positions values, respectively: where φ ij is the field value in the centre of cell j in the interpolation stencil of the vertex i, r ij is the centre of cell j in the interpolation stencil of the vertex i, and w ij is the weighting factor calculated as the inverse square distance between the position of vertex i and the centre of cell j. Finally, to obtain the unknown coefficient vector, C i , the weighted least-squares method is employed: where X is the n × 3 matrix whose row j is the position vector (r ij − r i0 ) of the cell j in the interpolation stencil of the vertex i, W is the diagonal matrix whose elements are the weighting factors for all cells in the interpolation stencil of the vertex i, and Φ i is the vector constituted by the elements (φ ij − φ i0 ) for all cells in the interpolation stencil of the vertex i.

Planar Extrudate Swell of Newtonian Fluids
The first benchmark case study that will be discussed is the planar extrudate swell of Newtonian fluids. Figure 3 shows a schematic representation of the computational flow domain, the boundary faces, and an indicative discretization mesh for the initial time-step (t = 0) and at steady-state. Cartesian coordinates are employed for the description of the planar flow domain, thus x = (x, y). The half width of the planar channel is denoted as H, which is considered as the scaling length. The inlet plane is taken sufficiently far upstream from the exit so that the flow is fully developed with a mean velocity U. Along the axis of symmetry the standard symmetry conditions are imposed. At the solid die wall the no-slip (tangential velocity is zero) and no penetration (normal velocity is zero) conditions are imposed. At the free-surface the kinematic, Equation (4), and dynamic conditions, Equations (6) and (7), are imposed. Finally, the outflow plane is taken sufficiently far downstream from the exit, so that the flow is uniform. The die exit of the planar domain is located at x = 5H from the inlet, and the outflow is located at x = 25H from the die exit.
In this section, we compare the results that were obtained with the newly-improved interface tracking algorithm with those given by Tanner [1], Georgiou and Boudouvis [27] and Mitsoulis et al. [28]. First, a mesh convergence sensitivity analysis is performed. Table 1 summarizes the main characteristics of the meshes employed in this preliminary study. The rectangular domain was initially discretized with 320 cells in the streamwise (x), 37 cells in wall normal (y), and one cell in span-wise directions, with this mesh being named M1. A linear stretch was employed in streamwise and wall normal directions to have the highest resolution at downstream edge of the die (at S) with the largest to smallest cell aspect ratio being equal to 2.5. An extensive investigation was carried out on the necessary grid spacing aspect ratio, and it was concluded that this value was efficient enough to track the free-surface. Five different meshes, designated M1, M2, M3, M4, and M5, were found to be sufficient to obtain accurate results for this problem. Table 1 shows the number of mesh elements used in each direction of the computational domain, the total number of mesh elements, and the degree of freedom (dof) of the numerical algorithm for each level of mesh refinement. Notice that each mesh is obtained from the preceding one by multiplying the number of cells in each direction by a factor of 1.5.   First, we compare our mesh convergence predictions of the steady-state swell ratio, χ, with the results found in the scientific literature [1,27,28]. The swell ratio is defined as the height of the free-surface away from the die exit, where the plug flow has been established, divided by the die height, i.e., χ = h 0 /H. For the mesh convergence analysis, we employed the Newtonian fluid at creeping flow conditions (Re = 0.1). Table 1 shows the extrudate swell ratio for the different mesh resolutions. The accuracy of the developed code is estimated via the application of Richardson's extrapolation [47] to the limit, by using the three most refined levels of mesh discretization. The extrapolated value of the swell ratio is equal to χ = 1.187, similar to the reference data of Mitsoulis et al. [28], which reported a value of χ = 1.186 using a FEM numerical algorithm, to the estimated extrapolated value of Tanner [1] with χ = 1.190 ± 0.002, and to the converged results obtained by Georgiou and Boudouvis [27] of χ = 1.186, as shown in Table 1. Additionally, notice that the relative error calculated between the extrudate swell ratio obtained using M5 (the most refined mesh) and the one extrapolated using the Richardson's technique [47] is equal to χ M5 − χ Extrapolated /χ Extrapolated × 100 ≈ 0.7%. Additionally, notice that, even when using M1 (the coarsest mesh), the relative error is only 2.3%. In terms of computational wall time, the simulation of the planar extrudate swell with M1 and M5 took approximately 0.94 and 24.3 hours, respectively, where all of the computations were performed in parallel using 80 processors on a computer with a 2.70 GHz Intel Xeon CPU E5-2680.
Subsequently, calculations for studying the robustness of the newly-improved interface tracking algorithm were pursued, by increasing the Reynolds number (Re ≤ 10) in the planar extrudate swell domain geometry. Figure 4 depicts the extrudate swell ratio, χ, as a function of the Re number, where the blue square symbols are the results obtained by the newly-improved interface tracking code using M5, and the dashed lines and red and green square symbols are the results found in the scientific literature [1,27,28]. As can be seen, the numerical results obtained are in very good agreement with the reference data.  [28], Georgiou and Boudouvis [27] and Tanner [1], respectively, and the blue square symbols represent the results obtained by the newly-improved interface tracking code. Figure 5 shows the transient evolution of χ for the planar extrudate swell domain geometry at Re = {0.1, 1, 2, 5, 7, 10}. The numerical results obtained for Re < 7 show an undershoot in the values of χ, before reaching the steady-state value. For Re ≥ 7, after the undershoot in the values of χ, an overshoot is present, and, ultimately, the steady-state value is reached, being approximately 0.99. This increase of inertia allowed for verifying the robustness of the improved interface-tracking algorithm, namely the Laplacian solver with least-squares volume-to-point interpolation, to handle abrupt changes in the mesh motion. The steady-state results for the primary field variables, magnitude of the velocity vector |u|, and pressure field p, in the form of contours, are shown in Figures 6 and 7, respectively, for the planar extrudate swell domain geometry, and at different Re numbers. From a detailed inspection of Figures 6 and 7, we can see that inertia substantially reduces the free-surface height, as already shown in Figure 4. This fact occurs, because inertia pushes the material to the center of the domain, generating a peak of negative pressure in the upper point near the die exit. This result was expected to occur due to the Poiseuille flow that developed in the upstream region of the domain geometry [48]. Additionally, for Re = 10, the contour of the magnitude of the velocity vector changes its behavior, when compared to the lower Re number cases. It can be seen that the maximum values of the velocity vector magnitude are extended after the die exit for the higher Re case. Finally, regarding the pressure contours, we can see that the increase of inertia does not change abruptly the pressure contours, and, at Re ≥ 7, we see an extension of the minimum pressure values from the top corner of the die exit to the center of the channel, which seems to cause the reduction in the free-surface height.  Finally, the efficiency of the newly-improved interface-tracking solver was assessed by comparing the required CPU wall time per time-step when using the original PISO algorithm or the enhanced consistent-PISO algorithm. Table 2 shows a comparison of the dimensionless time-step ∆t/(H 2 /ν) employed in the simulations of the planar extrudate swell and the elapsed time per time-step that is required by both PISO and consistent PISO algorithms for Re = {0.1, 1, 10}, while using mesh M5. For Re = 0.1, the consistent-PISO algorithm allows using a time-step 21.4 times higher than the one that is used by the PISO algorithm, and the CPU wall time per time-step that elapsed using the consistent-PISO algorithm is approximately half of the one taken by the PISO algorithm. For Re = 1 and Re = 10, the scenario is also favorable to the consistent-PISO algorithm to the detriment of the PISO algorithm, which allows for concluding that the enhanced correction procedure for the p − U fields is advantageous for efficiently simulating the extrudate swell problem. Table 2 also shows the total calculation time in hours when using M5 and the consistent PISO algorithm. As can be seen, all of the calculations finished in less than 1.5 days. Table 2. Comparison of the dimensionless time-step ∆t/(H 2 /ν) employed in the simulations and the CPU wall time (s) per time-step required by the PISO and consistent PISO algorithms for all Re, using mesh M5, for the planar extrudate swell flow of Newtonian fluids.

Axisymmetric Extrudate Swell of Newtonian Fluids
The second benchmark case study that is discussed in this work is the axisymmetric extrudate swell of Newtonian fluids. Figure 8 shows a schematic representation of the computational flow domain, the boundary faces, and an indicative discretization mesh for the initial time-step (t = 0) and at steady-state. Polar coordinates are employed for the description of the axisymmetric flow domain, thus x = (r, z). The half width of the axisymmetric channel is denoted as R, which is considered to be the scaling length. The boundary conditions imposed in the boundary faces are similar to the ones presented for the planar extrudate swell case study, with the exception that, in the axisymmetric domain, the two lateral boundary sides are considered to be wedge patches (i.e., the cylinder is specified as a wedge of small angle, e.g., 5º and one cell thick running along the plane of symmetry, straddling one of the coordinate planes), while, in the planar case study, they were considered as empty patches (i.e., this condition applies on each patch whose plane is normal to the third dimension for which no solution is required); and, instead of symmetryPlane at the bottom of the planar case, the axis of symmetry is considered as empty patch. The die exit of the axisymmetric domain is located at z = 5R from the inlet, and the outflow is located at z = 25R from the die exit.  Figure 9 illustrates the extrudate swell ratio, χ, for different Re with mesh resolution M5 (see Table 1) using the axisymmetric domain geometry. Similarly to the planar domain geometry, generally, the extrudate swell ratio that is obtained by the newly-improved interface tracking solver using the axisymmetric domain is in good agreement with the reference data [1,27,28].  [28], Georgiou and Boudouvis [27] and Tanner [1], respectively, and the blue circle symbols represent the results obtained by the newly-improved interface tracking code. Figure 10 shows the transient evolution of χ for the axisymmetric extrudate swell domain geometry at Re = {0.1, 1, 10}. Again, the numerical results obtained for Re < 7 show an undershoot in the values of χ, before reaching the steady-state value. When comparing with the planar case (see Figure 5), the magnitude of the undershoot is smaller for the present case, which can be attributed to the round boundary surface of the axisymmetric domain geometry. For Re ≥ 7, an almost imperceptible undershoot is present when using the axisymmetric domain geometry (contrarily to the planar domain geometry), followed by an overshoot, and ultimately the steady-state χ value is reached, being approximately 0.97 (2% less than the value that is obtained for the planar domain geometry). The steady-state results for the primary field variables, magnitude of the velocity vector |u|, and pressure field p, in the form of contours are shown in Figures 11 and 12, respectively, for the axisymmetric extrudate swell domain geometry and at different Re numbers. From a detailed inspection of Figures 11 and 12, we can also see here that inertia substantially reduces the free-surface height, as already shown in Figure 9. Additionally, when compared to the lower Re number cases, we can see that, only for Re = 10, the contour of the magnitude of the velocity vector changes its behavior, where the maximum values are extended until the die exit for the latter case. Again, the increase of inertia, does not change abruptly the pressure contours, and only at Re = 10, we see an extension of the minimum pressure contours from the top corner of the die exit to the center of the channels, which seems to cause the reduction in the free-surface height (but with a lesser effect in the axisymmetric domain geometry).
When comparing the contour results from both planar and axisymmetric extrudate swell domain geometries, we can state that the maximum velocity vector magnitude is obtained for the axisymmetric domain, being two times the magnitude of the inflow average velocity, in contrast with the planar case, where it is only 1.5 times higher than the magnitude of the inflow average velocity. This result was expected to occur due to the Poiseuille flow, which is developed in the upstream channel of both domains [48]. Finally, regarding the pressure contours, we can see that, for the axisymmetric domain, the computed maximum and minimum pressure values are symmetric, as expected, due to domain symmetry, in contrast with the planar domain case, where the minimum pressure value that is obtained is 3.5 times lesser than the maximum one. In both planar and axisymmetric cases, as the Re increases the extrudate swell ratio and pressure decreases, which physically states that the inertia forces stretch out the material and prevent the swelling, with an immediate effect on the reduction of pressure losses [28]. Finally, Table 3 shows a comparison of the dimensionless time-step ∆t/(H 2 /ν) employed in the simulations of the axisymmetric extrudate swell and the CPU wall time (s) per time-step required by both PISO and consistent PISO algorithms for Re = {0.1, 1, 10} using mesh M5. For Re = 0.1, the consistent-PISO algorithm allows for using a time-step 21.0 times higher than the one used by the PISO algorithm, and the CPU wall time per time-step elapsed using the consistent-PISO algorithm is 0.58 times the one taken by the PISO algorithm. Again, for Re = 1 and Re = 10, the scenario is also favorable to the consistent-PISO algorithm in detriment of PISO algorithm. Additionally, notice that, when comparing with the values that are shown in Table 2, the ratios of CPU wall time per timestep (C-PISO/PISO) are lower for the planar extrudate swell, which mean that the CPU wall time gains using the C-PISO algorithm seems to be higher for non-smooth geometries. Table 2 also shows the total calculation time in hours when using M5 and the consistent PISO algorithm. As for the planar extrudate swell, all of the calculations finished in less than 1.5 days.

Conclusions
A numerical formulation for efficient moving mesh interface tracking simulations of free-surface flows was presented and implemented using the finite-volume method. The implementation was performed in the open-source OpenFOAM framework [37], where the interface is tracked in a semi-implicit manner inside the consistent second-order timeaccurate non-iterative Pressure-Implicit with Splitting of Operators (consistent PISO) algorithm for the numerical solution of incompressible fluid flows. Additionally, the moving mesh was adjusted to the time varying shape of the interface, using a Laplacian scheme with a least-squares volume-to-point interpolation, which allowed for robust and stable deformations of the interface.
The improved algorithm was assessed in terms of the accuracy and efficiency for the fluid flow simulations in planar and axisymmetric extrudate swell of Newtonian fluids. A mesh sensitivity analysis allowed for obtaining a grid refinement level from which the calculated extrudate swell ratio differs from the extrapolated value of 0.7%. Subsequently, the robustness of the numerical algorithm was pursued, increasing the Reynolds number from Re = 0.1 to Re = 10. The extrudate swell ratio that was obtained for both domains compared well with the results found in the scientific literature for that range of Re. Additionally, the contours for the magnitude of the velocity vector and pressure fields are also shown, and a detailed study of the contours reveals that the obtained results are physically meaningful. Finally, the efficiency of the improved numerical solver was evaluated by comparing the CPU wall time (s) per time-step for both the PISO and consistent PISO algorithms. The results obtained allowed for concluding that the consistent-PISO is, at maximum, 47% and 42% faster than the PISO algorithm for the planar and axisymmetric extrudate swell domain geometries, respectively.
In summary, the results presented here show that the newly-improved interface tracking code, developed using an open-source framework, can accurately and efficiently predict the Newtonian extrudate swell. The code that was implemented here is being currently extended to handle viscoelastic fluid flow calculations and non-isothermal processes. For the viscoelastic fluid flow calculations, we will use the quasi-linear Oldroyd-B and exponential PTT rheological models. The former will be used due to the numerical instabilities that are caused by the infinite polymeric stresses generated at singular points, which will verify the robustness of the numerical implementation; and the latter will be used because it is more suitable for approximating the behavior of polymer melts, where the extensibility parameter introduces elongational and shear thinning in the fluid model. Additionally, the solvent and polymeric viscosities and the relaxation time of the viscoelastic fluid will be considered as temperature dependent, by employing the Williams-Landel-Ferry relation.