Spectral/hp element methods for plane Newtonian extrudate swell

Spectral/hp element methods and an arbitrary Lagrangian-Eulerian (ALE) moving-boundary technique are used to investigate planar Newtonian extrudate swell. Newtonian extrudate swell arises when viscous liquids exit long die slits. The problem is characterised by a stress singularity at the end of the slit which is inherently difficult to capture and strongly influences the predicted swelling of the fluid. The impact of inertia (0<Re<100) and slip along the die wall on the free surface profile and the velocity and pressure values in the domain and around the singularity are investigated. The high order method is shown to provide high resolution of the steep pressure profile at the singularity. The swelling ratio and exit pressure loss are compared with existing results in the literature and the ability of high-order methods to capture these values using significantly fewer degrees of freedom is demonstrated.


Introduction
In this article, we investigate the extrudate swell phenomenon, which is a radial swelling of free liquid jets exhibited by viscous fluids exiting long die slits. This jet swelling is particularly strong for viscoelastic fluids but is also exhibited by low Reynolds number Newtonian fluids. The prediction of the swelling ratio is very important in a range of industrial processes such as inkjet printing, extrusion moulding or cable coating.
The swelling of Newtonian jets is mainly characterised by the reorganisation of the velocity profile from the parabolic Poiseuille flow inside the die to plug flow downstream (Tanner [2002]). This transition is characterised by the sudden jump in the shear stress at the die exit (Russo [2009]). Inside the die, the shear stress at the wall is at its maximum with particles sticking to the wall (for the no-slip boundary condition). Then immediately after the die exit, the removal of the wall shear stress causes a boundary layer to form at the free surface. In this layer, the parabolic velocity profile adjusts itself so as to satisfy the condition of zero shear stress at the free surface. This sudden jump in the shear stress at the die exit causes an almost instantaneous acceleration of the particles at the free surface causing the fluid jet to swell.
Due to the presence of this stress singularity at the die exit, numerical simulations of the extrudate swell phenomenon are particularly challenging. Analytically, this singularity originates from the sudden change in the boundary condition from the wall of the die to the free surface of the exiting jet. This "jump" in the boundary condition yields steep and infinite stress and pressure concentrations at the singular point. These infinite stress values near the singularity affect the accuracy of the numerical solution and the size of the swelling and therefore need to 1 arXiv:1408.5167v1 [math.NA] 21 Aug 2014 be resolved as accurately as possible. In this contribution, we use a spectral/hp element method to improve our ability to capture these stress concentrations. Traditional discretisation methods such as finite differences or low-order finite elements require a very large number of degrees of freedom to resolve these sharp stress variations.
In this article, we will describe a spectral/hp method that is capable of approximating the infinite stress values with an exponential increase in the extreme values of the pressure with p-refinement. This demonstrates that our high order method provides a high-quality approximation of the stress singularity with a very low number of degrees of freedom. We will give detailed information about the pressure and velocity in the vicinity of the singularity for a wide range of Reynolds numbers (0 ≤ Re ≤ 100) and for slip along the die wall. We demonstrate that our method predicts swell ratios and exit pressure loss corrections in excellent agreement with a recent numerical study of Mitsoulis et al. [2012] for our coarsest approximation P = 10. Mitsoulis et al. [2012] used a low order finite element method with a high mesh refinement around the singularity.
Typically, a decrease in the swelling is observed for an increase in the resolution of the singularity. In the existing literature, the stress values at the singularity are rarely addressed. Salamon et al. [1995] investigated the role of surface tension and slip on the singularity numerically and analytically. They demonstrated that a very fine mesh near the singularity is needed to predict the singular pressure and stress behaviour with sufficient accuracy. Georgiou and Boudouvis [1999] compared the singular finite element method with the regular finite element method for the extrudate swell problem. In the singular finite element method basis functions in the elements around the singularity are enriched with the local asymptotic solution for the singularity. They demonstrated that with this method the predictions of the swell ratio converged. However, the singular finite element method requires the correct asymptotic behaviour of the pressure at the corner singularity and the asymptotic solution for the pressure is obtained assuming Stokes-like behaviour around the singularity. This means this approach is only accurate for Re = 0. Indeed, Georgiou and Boudouvis [1999] found that the singular finite element method was outperformed by the regular finite element method for extrudate swell including inertia. Our method is capable of resolving the stress singularity with spectral convergence properties without making any assumptions on the form of the singularity.
Inertialess extrudate planar Newtonian swell has been investigated in terms of swell ratios using low order finite elements by a wide range of authors (Tanner [1973], Nickell et al. [1974], Crochet and Keunings [1982], Reddy and Tanner [1978]). Tanner [2002] provides a review of inertialess Newtonian swell ratio results. Only very few investigations involved the use of higher order methods. Ho and Rønquist [1994] provided the first extrudate swell computation with a spectral method for one coarse mesh with 8 spectral elements with polynomial order 4 for Re = 0. They predicted a swell ratio of 1.1840. Russo [2009] used the spectral element method to predict free surface profiles and swell ratios for 0 ≤ Re ≤ 10 and surface tension for 4 spectral elements with polynomial order 6 < P < 14. We will use a spectral element mesh with 14 spectral elements and 10 ≤ P ≤ 16 with a smaller element size around the singularity providing a much higher resolution there compared with previous studies. We will provide results for 0 ≤ Re ≤ 100 and for a slip condition along the die wall.
The paper is organised as follows. In Section 2, we will introduce the governing equations for the description of Newtonian free surface flow and the equations of motion for the mesh movement. We will conclude this Section with a description of the boundary conditions for the extrudate swell problem and the definition of the quantities of interest such as swelling ratio and exit pressure correction. In Section 3, we describe the numerical discretisation of the governing equations. In Section 4, we give numerical results for the impact of inertia and slip on the extrudate swell problem including detailed plots for velocity and pressure profiles in different parts of the domain.

Governing Equations of the Fluid
The free surface motion of an incompressible fluid flow can be characterised by the incompressible Navier-Stokes equations describing the motion of the fluid and the motion of the free surface. On a moving domain Ω t ⊂ R d , t ∈ I ≡ (t 0 , T ), they can be expressed as where u is the velocity, p is the pressure, D = 1 2 (∇u + ∇u T ) is the rate of deformation tensor, Re is the Reynolds number, u 0 is the velocity field at t = t 0 and u D is the assigned Dirichlet boundary condition.
The motion of the free surface, Γ f , is characterised by the following boundary conditions where w is the velocity of the free surface, σ is the surface tension coefficient, κ is the curvature of the free surface, n is the unit outward normal to the free surface and [σ σ σ] denotes the jump in the Cauchy stress tensor across the free surface. In order to track the free surface motion computationally, the grid points of our computational mesh at the free surface are moved with the normal fluid velocity, which ensures that particles do not cross the interface and therefore that the kinematic condition (2a) is satisfied.
To avoid mesh distortion, the mesh points in the interior of the domain are moved with an arbitrary speed. This use of arbitrary mesh movement is known as the arbitrary Lagrangian-Eulerian (ALE) technique. The ALE formulation relates the Navier-Stokes equations on the moving domain (1) to a formulation on a referential configurationΩ t 0 . At each t ∈ I, each point of the reference configuration Y is then associated to a point x x x in the current domain Ω t using the so-called ALE-map [Donea et al., 2004, Scovazzi and Hughes, 2007, Pena, 2009, Nobile, 2001, that is, where Y is called the ALE coordinate and x x x is the Eulerian coordinate. The movement of the mesh, can then be characterised by the following quantities 1. the mesh velocity 2. the material time derivative in terms of the time derivative with respect to the ALE-frame Equation (1) in the ALE-formulation reads is the rate of deformation tensor in the Eulerian frame of reference.

Governing Equations of the Mesh
In addition to the motion of the fluid, we need to find a sensible way to describe the domain movement. In general, the domain movement is characterised by the movement of its boundary ∂Ω t and can be described using the domain or mesh velocity w (Ho and Rønquist [1994], Robertson et al. [2004]), the ALE-mapping R(t) (Nobile [2001], Pena [2009]) or the displacement d = ∆tw (Choi and Hulsen [2011]). In the present work, we describe the domain movement using the mesh velocity, w. For the domain movement, we choose boundary conditions such that the kinematic boundary condition is satisfied and mesh distortions are kept to a minimum, that is, where s is the unit tangential vector on the free surface boundary. In order to guarantee smooth mesh movement in the interior, we solve an elliptic problem for the mesh velocity, given by ∆w = 0 on Ω(t).
subject to the boundary conditions (7). This harmonic mesh movement preserves a high quality mesh for small displacements and has been employed, for instance, by Ho and Rønquist [1994], Nobile [2001] and Pena [2009]. However, for higher mesh deformations, other elliptic problems may be solved for the movement of the domain, such as elliptic operators arising from Stokes or elasticity problems (see the monograph of Deville et al. [2002] for further details).
Inflow: Figure 1: Schematic of the symmetric die swell flow configuration. L 1 is the length of the die which has fixed boundaries and is of half-height H. L 2 is the length of the outflow region, the boundaries of which are free to move. Boundary conditions are provided for each surface.

Computational Domain and Quantities of Interest
Consider the extrusion of a Newtonian liquid from a planar die. The schematic diagram of the employed planar die geometry is depicted in Figure 1. We consider a die of length L 1 and height H, and an exit region of length L 2 . The length of the die is chosen sufficiently long in order to guarantee a fully developed flow far upstream of the exit plane. In the following, we pay special attention to the following two quantities of interest: the swelling ratio and the pressure exit correction factor. In practice, the extrudate swell ratio is of importance in extrusion processes and the excess pressure loss gives an indication as to how much extra pressure has to be applied to achieve certain swell ratios. The swelling ratio, χ R , is defined as where H is the half-height of the die and h f is the half-height of the liquid jet at the outflow boundary. The swelling ratio is a function of several parameters where u is the average inflow velocity, Re is the Reynolds number and B sl is the slip parameter along the die wall. The dimensionless pressure exit correction factor, n ex , is defined as where ∆p is the pressure drop between the inlet and the outlet plane, ∆p 0 is the pressure drop between the inlet and the exit of the die for fully developed Poiseuille flow and σ w is the shear stress at the channel wall corresponding to fully developed Poiseuille flow. Here, the pressure differences are taken along the centreline. In particular, the pressure drops are given by (Tanner [2002]) In our computations, we employ the following boundary conditions as depicted in Figure 1 for a half-channel height of H = 1. We assume the flow is symmetric and along the symmetry line, we set v = 0 and σ xy = 0. Note that, σ xy = 0 is set through the boundary integral in the momentum equation (25). For the die swell geometry this means that there is no contribution of the Neumann boundary integral in the momentum equation along the symmetry line. At the die wall, we either impose no-slip boundary conditions (i.e. u = 0) or Navier's slip condition. The latter is a mixed boundary condition of Dirichlet and Neumann type. For the extrudate swell geometry depicted in Figure 1, we set v = 0 and impose σ xy = 1 B sl u through the Neumann boundary term in the momentum equation (25). This means for the velocity component u along the slip boundary Γ sl , we obtain the boundary integral where e x is the unit vector in x-direction. At the outflow, we employ an open outflow boundary condition. We assume a reference pressure p = 0 along the outflow boundary and the remaining terms in the Neumann boundary integral along the outflow boundary in the momentum equation are evaluated along with the volume integrals. At inflow, we either impose the parabolic profile in combination with no-slip along the die wall or the profile (Kountouriotis et al. [2013]) in combination with the slip boundary condition. In the extrudate swell problem the velocity field undergoes a transition from Poiseuille flow inside the die to plug flow in the free jet. Due to the conservation of energy the flow rate in the die has to be the same as in the uniform plug flow, which yields where h plug is the height of the fluid jet in the uniform flow region and u(y) is the parabolic Poiseuille flow profile. We have 0 < u plug < u max , which means that while particles at the free surface accelerate when exiting the die the flow near the centreline decelerates.
For the mesh velocity, we employ the following boundary conditions. We consider the mesh to be fixed at inflow, the die wall and along the symmetry line, i.e. w = (w x , w y ) = 0. At the outflow boundary, we allow the mesh to move in the y-direction, i.e. ∇w y · n = 0, and fix it in the x-direction, w x = 0. At the free surface, we enforce the kinematic boundary condition through the mesh velocity in terms of a Dirichlet boundary condition for the mesh-velocity, i.e.
To avoid mesh distortion, we choose to move the mesh along the free surface boundary only in the y-direction. The mesh is moved with sufficient velocity w y into the y-direction to ensure that no particle crosses the interface, that is, 6 3 Numerical Discretisation

Spectral Element Discretisation
Consider the decomposition of the domain Ω δ t into N el non-overlapping elements. These elements are each mapped to a standard element on which the unknowns are approximated using a modal polynomial expansion basis proposed by Dubiner [1991] and extended by Karniadakis and Sherwin [2005] given by Here, φ 0 and φ P are the linear finite element basis functions and is the usual quadratic hierarchical expansion mode for quadratic elements. Furthermore, P denotes the highest polynomial order of the hierarchical expansion and P (α,β) p (ξ) denotes the p th -order Jacobi polynomial.
Two-dimensional functions u(x, t) can be approximated on two-dimensional standard quadrilaterals, defined as Ω st = {−1 ≤ ξ 1 , ξ 2 ≤ 1}, using a tensor product of the one-dimensional modal expansion basis functions φ p , that is, with the reference coordinates given by involving the inverse of the mapping χ χ χ e . Here, the mapping, χ χ χ e , between the local coordinates (ξ 1 , ξ 2 ) and the physical coordinates (x, y) approximates the geometry with the same order polynomial space as the solution, that is, Details on the construction of this mapping can be found in Karniadakis and Sherwin [2005]. 7

Weak Formulation
Introducing the function spaces in the current frame with respect to the reference configuration the weak formulation of the system of equations (6) leads to the following problem definition.
Problem 3.1 (Weak formulation of incompressible Navier-Stokes equations). For almost every where Γ N (t) is the Neumann boundary and Γ f (t) is the free surface boundary.
We choose the same trial and test function space for the mesh velocity as for the fluid velocity, i.e. we choose and we solve Equation (8) with the boundary conditions (7) using a continuous Galerkin method. The weak formulation for the mesh movement can therefore be expressed as below.
subject to the boundary conditions (7).
The position of the new nodes of the mesh can be obtained via Equation (4), that is, 8

Discrete ALE formulation
As mentioned above, we have two referential domains to consider in the ALE formulation. Firstly, let Ω δ t be the union of all non-overlapping mesh elements in the Eulerian frame at time t and secondly, letΩ δ t 0 denote the union of all mesh elements in the referential frame. Consider the following discrete trial and test function spaces for the fluid and mesh velocities and for the pressure field. Alternatively, these spaces can be expressed as (see Pena [2009]) Here, P c P (Ω δ t 0 ) denotes the globally continuous space of polynomials of degree P over the reference mesh, that is, P P (Ω δ t 0 ) denotes the space of piecewise continuous polynomials of degree P over the reference mesh, that is, P c P (Ω δ t ) denotes the globally continuous polynomial space over the Eulerian mesh and P P (Ω δ t ) denotes the piecewise continuous polynomial space over the Eulerian mesh. Here, g δ Ω e t 0 denotes the restriction of g δ to the spectral element Ω e t 0 , P P (Ω st ) is the space of polynomials of degree P defined on the standard element given by the expansion basis (20). Note that, the pressure is discretised with polynomials of order 2 lower than the velocities to satisfy the LBB condition [Brezzi, 1974]. The spaces V δ D (Ω δ t ) and Q δ (Ω δ t ) include the discrete ALE mapping, which can be expressed as [Nobile, 2001] involving the geometrical mappings, χ χ χ e (t), at each t, from the standard element Ω st to each element Ω e t , that is, 9 wherex pq (t) denotes the expansion coefficients at time t and the iso-parametric mapping, χ χ χ e (t 0 ), from Ω st to Ω e t 0 , defined as Using these space definitions and an implicit Euler time-integration scheme, the semi-discrete Navier-Stokes equations are expressed as follows.
Problem 3.3 (Semi-discrete Navier-Stokes ALE formulation). For each n, let t n = t 0 + n∆t, ). Here, we linearise the convective term in the momentum equation by setting u * δ = u n δ , which is an extrapolation of the velocity of the same order as the implicit Euler scheme. Note that, the index S for normals and curvature in the boundary integral over Γ f (t n+1 ) indicates that these quantities are determined from a cubic spline representation of the free surface according to equation (48) and (49) defined below.

Matrix formulation
The discrete ALE formulation involves the following matrices and a modified Helmholtz matrix The equation system (39)-(40) can then be written for each element in algebraic form as whereû g andp g are the vectors of unknown global coefficients, H g , D g = (D x 1 , D x 2 ) are the global matrices assembled from the elemental matrix contributions. The resulting system of equations is then solved using a multi-level static condensation technique introduced by Ainsworth and Sherwin [1999], Sherwin and Ainsworth [2000] and Karniadakis and Sherwin [2005] for the Stokes equations in fixed domains.

Discretisation of Mesh Movement
Even though solving Problem 3.2 yields continuous mesh movement, the free surface boundary might not be sufficiently smooth. The free surface boundary undergoes the largest deformation and its movement involves the evaluation of outward normals, n, in Equation (7), across multiple elements. Note that, a standard Galerkin method with a C 0 -continuity across elements is not sufficient to determine a well-defined normal at element edges. To alleviate this problem, we represent the free surface using a cubic spline, S(x, t) ∈ C 2 (Γ f ) to ensure sufficient smoothness of free surface boundary edges of the mesh. The cubic spline can then be used to determine the unit outward normals n and the curvature κ of the free surface using These expressions are then used to evaluate the free surface boundary condition for the mesh velocity given by Equation (7) and the free surface boundary integral in the momentum equation For given u n , we perform the mesh movement in the following way. First, we determine the cubic spline through all the quadrature points along the free surface. Let (x i , y i ), 1 ≤ i ≤ N, be the physical coordinates of the N quadrature points along the free surface. Then, we construct a cubic spline S(x, t) = S i (x, t) for each where we enforce continuity and smoothness We employ the not-a-knot boundary condition on the spline, that is, We then solve the elliptic problem (28) using the continuous Galerkin method, determining where L g is the global Laplace matrix given by subject to the boundary conditions (7), which include the normal determined by the cubic spline according to (48). The mesh velocity resulting from the solution of Equation (56), denoted byw, is then used to update the coordinates of the mesh nodes using a third order Adams-Bashforth-Scheme for Equation (4), that is, This equation is solved pointwise in the strong form for each quadrature point. However, in practice, we do not move all the mesh nodes of every element. We only move all the quadrature points along the free surface boundary introducing curved edges along the free surface boundary.
In the interior of the domain, we just move the corner vertices of every element keeping the interior edges of the domain straight.
Using the new coordinates of all mesh nodes, we compute the mesh velocity at the new time level pointwise as 12

Mesh Configuration
We use a mesh consisting of N el = 14 elements as shown in Figure 2 and refine the mesh by increasing the polynomial order P . We consider a die of length L 1 = 10 and an exit region of length L 2 = 10. The entry length is sufficient to guarantee a fully developed flow far upstream from the exit of the die. The exit length is chosen sufficiently long to allow the free surface to reach a constant downstream height for a large range of Reynolds numbers. For high Reynolds numbers, the free jet length might be insufficient to guarantee a fully developed plug flow profile at outflow. However, the use of the open outflow boundary condition enables us to predict the correct swelling ratios truncated at the outflow boundary location (see Mitsoulis et al. [2012]). Throughout this section, we choose a time step of 5 × 10 −3 .

Numerical results for Re = 0
Inertialess Newtonian extrudate swell has been investigated in a number of publications. Table 1 summarises some of the swelling ratios obtained by a range of authors for plane Newtonian die swell. Tanner [2002] used the results in the literature to estimate an extrapolated value for planar die swell of χ R = 1.190 ± 0.002. In general, an increase in the number of degrees of freedom yields less swelling.   Table 2 lists a comparison of the pressure exit correction for Re = 0 of our scheme and the swell ratio for increasing mesh refinement with the results obtained by Taliadorou et al. [2007]. We obtain close agreement for a much smaller number of degrees of freedom, which demonstrates that p-refinement is effective for the Newtonian extrudate swell even though the result is polluted by Gibbs oscillations in the pressure around the singularity (Figure 3c). The Gibbs oscillations in the pressure stay confined to the elements adjacent to the singularity. Increasing the Reynolds number leads to a dampening in the oscillations in the elements adjacent to the singularity and the extreme values of the pressure at the singularity decrease significantly (Figure 11b). As shown in Figure 3c increasing the polynomial order yields an increase in the number of oscillations. However, the amplitude of each oscillation is reduced with increasing polynomial order P . Increasing the polynomial order also has the effect of exponentially increasing the maximum value of the pressure and sharply increasing the minimum value of the pressure at the singularity which reflects an improved approximation of the infinite pressure value at the singularity (Figure 3d). While the infinite pressure values at the singularity hamper the rate of  convergence of the numerical pressure solution, the values of the velocity components along the free surface are converged for P ≥ 10 (see Figure 3(a), (b)).

Impact of inertia
Inertia causes a decrease of the swelling and the liquid jet eventually contracts for sufficiently high Reynolds numbers. We performed computations for Reynolds numbers ranging from 0 to 100. We start by computing the extrudate swell for Re = 0 and initialise this computation with the solution of the corresponding stick-slip problem. After having obtained the extrudate swell for Re = 0, we increase the Reynolds number in steps of 1 from 1 to 10 and in steps of 10 from 10 to 100, each time using the result of the converged extrudate swell of the previous lower Reynolds number as the initial condition. As the convergence criterion, we choose a change of the maximum absolute value of all variables including the mesh velocity of less than 10 −6 . Figure 4 and Table 3 shows the comparison of the swelling ratios obtained with our algorithm with the results of Mitsoulis et al. [2012], which are in excellent agreement. Figure 5 displays the corresponding free surface spline profiles. We observe that the swelling ratio decreases at an accelerating pace with increasing Reynolds number until Re = 6. For Re = 6, we see the onset of a delayed die swell in which the fluid surface first goes through a minimum before it swells again. The delay in the swelling of the jet increases with increasing Reynolds number from Re = 6 to Re = 10. For Re = 9 and Re = 10, the fluid contracts (χ R < 1) but still experiences some swelling after going through a minimum near the die exit. For Re = 20 to Re = 100 the fluid does not experience any delayed swelling and contracts. For 10 < Re < 40 the fluid contracts very fast with increasing Reynolds number. This trend in the contraction rate with increasing Reynolds number then slows down and approaches a limit for 40 < Re < 100. The limit for infinite Reynolds number was estimated by Tillett [1968] who performed a boundary layer analysis for a free Newtonian jet and predicted a limiting value of χ R = 0.8333 for infinite Reynolds number.
We explore the contour plots of the velocity field for a range of Reynolds numbers in Figures 6 (horizontal velocity component u), 8 (vertical velocity component v). With increasing Reynolds number the horizontal velocity increases along the centreline, the vertical velocity near the singularity induced by the sudden change in the boundary condition decreases and the transition zone under the free surface from Poiseuille flow in the die to plug flow is extended downstream. This shows that with increasing Reynolds number the particles along the centreline are accelerated and decelerated near the free surface yielding the contraction of the free fluid jet. This is indeed the behaviour we would expect as particles leaving the die will deviate less from their initial path for increasing inertia. As pointed out by Mitsoulis et al. [2012] in order to accommodate the whole transition zone the domain length of the free fluid jet should be chosen as L 2 = Re . However, we employ open boundary conditions at outflow which enable us to compute the extrudate swell accurately in the truncated domain with L 2 = 10. As demonstrated by Mitsoulis and Malamataris [2011] the results for extrudate swell with a domain length L 2 = 6 are virtually identical with those from long domains with L 2 = Re , for all variables, when using the open boundary condition at outflow. However, in this case, the swell ratio results are only correct up to the truncated length as they continuously drop beyond the truncated domain. A small discrepancy between swell ratios for different domain lengths can therefore be expected.
To investigate the transition from Poiseuille flow to plug flow for increasing Reynolds number further, we plot the velocity and pressure along different paths in the domain. Figure 7 displays the velocity components along the symmetry line (i.e. v = 0) and along the free surface boundary. In Figure 7a, we see the smooth transition of the velocity field from the maximum of the parabolic profile to the average plug flow velocity given by Equation (17), i.e. u plug = 1/χ R . As the swell decreases with increasing Reynolds number the plug flow value of the velocity increases with increasing Reynolds number. With increasing Reynolds number the change from the maximum parabolic value of the velocity component u to the plug flow value shifts further downstream. For Re = 0, the velocity reaches the plug flow value at around x ≈ 3, for Re = 10 at x ≈ 6 and for Re = 50 the plug flow value is not reached within our computational domain. However, as pointed out above, due to the use of open boundary conditions at outflow, the velocity and pressure profiles stay accurate even if they are truncated at outflow.
Along the free surface boundary (Figure 7b), the velocity component u increases sharply near the die exit until it reaches the plug flow value while the velocity component v goes through a maximum near the die exit for Re = 0 and Re = 3 and through a minimum for Re > 7, when particles are no longer constrained by the no-slip boundary condition (Figure 7c). This causes the swell (for v > 0) or the contraction (for v < 0) of the free surface near the die exit until the surface is sufficiently curved to obtain a zero total shear stress (i.e. t · σ σ σ · n = 0). Further downstream when the free surface boundary has reached its maximum swelling value, the vertical velocity component reaches zero in accordance with the condition of no particle penetration along the surface (horizontal free surface boundary has outward normal n = (0, 1) and therefore u · n = v = 0). The maximum value of v along the free surface decreases with increasing Reynolds number (0 ≤ Re ≤ 5). For the range of Reynolds number that causes a delayed die swell the velocity component v first undergoes a sharp minimum and then goes through a maximum (6 ≤ Re ≤ 10). For the range of Reynolds numbers that cause a contraction of the free Newtonian jet, the velocity component v goes through a minimum and then slowly approaches zero (Re > 10). Figure 9 shows the velocity components in the cross stream wise direction at inflow (x = −10), near the die exit (x = −0.2, x = 0.2), further downstream in the free jet region x = 1 and at outflow x = 10. The velocity component u, is parabolic at inflow, shortly before the die exit (x = −0.2) the parabolic profile flattens inside the die, after the die exit the parabolic profile flattens further and builds a boundary layer in which it goes through a minimum x = 0.2, then flattens increasingly until the plug flow value is reached. The vertical velocity component, which is zero at inflow, forms a parabolic-like profile with a small boundary layer near the die exit inside the die, which first sharpens shortly after exiting the die and then relaxes back to the zero value.
In the contour plots for the pressure p displayed in Figure 10, we observe that the pressure isobars are curved near the die exit and in the free jet region into the downstream direction for low Reynolds number (Re = 0, 3, 7) and into the upstream direction for higher Reynolds numbers (Re > 10). The change in the pressure becomes more apparent when we explore the pressure values along the symmetry line (Figure 11a). Inside the die, the pressure gradient is constant as expected for Poiseuille flow. However, near the die exit (x = 0) the pressure smoothly approaches zero for the plug flow. For higher Reynolds numbers the pressure on the centreline goes through a minimum. This behaviour of the pressure yields a shift in the pressure values at inflow, which is expressed by the pressure exit correction as defined in Equation (11).

Impact of slip
To alleviate the pressure singularity at the die exit, we investigate the effect of slip along the die wall on the dependent variables for Re = 0. We therefore change the inflow profile according to Equation (16) and employ the slip condition (14) along the die wall. We explore the velocity field and the pressure along the free surface for slip parameter values of B sl = 0.01, B sl = 0.1 and B sl = 0 (no-slip) in Figure 12. With the introduction of slip along the wall, the horizontal velocity component experiences a smooth transition at the die exit in vast contrast to the kink at the singularity that is observed for the no-slip condition (B sl = 0) along the wall (Figure 12a). The change for the vertical velocity remains sudden and features a kink at the singularity. However, the maximum value of the vertical velocity component decreases with increasing slip (Figure 12c).
The pressure profile at the singularity is changed drastically with slip along the wall and the Gibbs oscillations disappear (Figure 12e, (f)). Even though the minimum of the pressure does not show a converging trend in the range of the employed polynomial orders, its value only increases slightly with increasing P (Figure 12f). Table 4 lists the swelling ratios for increasing polynomial order, P , for B sl = 0.1 and B sl = 0.01. The swelling ratios are converged to three decimal places. Figure 12b and 12d show that the velocity values are converged for P ≥ 10. The free surface spline for increasing slip parameter is shown in Figure 13. Increasing the slip parameter yields a decrease in swelling.

Conclusions
In this article, we have demonstrated the capabilities of the high-order spectral element method in the resolution of the stress singularity at the die exit in the plane Newtonian extrudate swell problem. We have shown that the spectral method approximates the infinite pressure value with exponentially increasing extreme values for increasing polynomial order. This high resolution approximation of the steep stress profiles yields excellent predictions of the swelling ratio. Our method predicts the same swelling ratio in comparison to low order finite element methods with significantly fewer number of degrees of freedom. The only drawback of our high order method is the Gibbs oscillations, which appear in the vicinity of the singularity for the pressure approximation. These Gibbs oscillations are intrinsic to high order methods and they occur in the approximation of discontinuities or steep profiles. However, we have demonstrated that for the extrudate swell problem, the Gibbs oscillations stay confined to one element next to the singularity and their amplitude decreases significantly with increasing polynomial order. This small pollution in the pressure profile is the price to pay in the high order method for the otherwise excellent prediction of the steep pressure increase at the singularity.
We have given detailed results for a wide range of Reynolds numbers 0 ≤ Re ≤ 100 in terms of swell ratios, exit pressure losses, free surface profiles and velocity and pressure values. For the free surface profiles, we find three extrudate swell regimes. The first is a reduction in swelling (Re ≤ 6), the second is a regime of a delayed swelling (7 ≤ Re ≤ 10) and the third a contraction of the free liquid jet (10 < Re ≤ 100). With increasing Reynolds number the maximum pressure values decrease and the Gibbs oscillations decrease. We have then investigated the effect of slip along the die wall. We have observed a reduction of the swelling for different slip parameters B sl = {0.01, 0.1} and have observed a drastic change in the pressure profile which showed no occurrence of Gibbs oscillations.