Transient 3D finite element method for predicting extrudate swell of domains containing sharp edges

A new transient 3D finite element method for predicting extrudate swell of domains containing sharp edges is proposed. Here, the sharp edge is maintained over a large distance in the extrudate by describing the corner lines as material lines. The positions of these lines can be used to describe the transverse swelling of the 2D free surfaces and expand the domain over which a 2D height function on the free surfaces is applied. Solving the 2D height functions gives the positions of the free surfaces. First a 2D axisymmetric case was tested for comparison, using three different constitutive models. The Giesekus, linear Phan-Thien Tanner (PTT) and exponential PTT constitutive models all showed convergence upon meshand time-step refinement. It was found that convergence remains challenging due to the singularity at the die exit. The new method is validated by comparing the final volume change of the extrudate of a 3D cylinder to the final volume change of a reference mesh of the 2D axisymmetric case. Finally, simulations were performed for different, complex, die shapes for a viscous fluid and for viscoelastic fluids. The results compared favorably with literature. Viscoelastic results, using the Giesekus model and the exponential PTT model, were compared for different Weissenberg numbers and different values for the non-linear parameters of the constitutive models. It was found that the swell is highly dependent on the rheological parameters and the constitutive model used.


Introduction
Extrusion is a widely used process to create products with a fixed cross-sectional profile. Many applications require cross-sections of complex shapes, containing sharp edges. Common requirement on the extrudate is dimensional precision. The dimensions of the extrudate, however, are highly influenced by a phenomenon called extrudate swelling, where the extrudate will start to expand due to internal stresses once it leaves the die. When the extrudate has such a shape that it contains sharp edges, the material lines exiting from the corners will swell in two directions. Also, the sharp kink in the surface will be maintained over a large distance from the exit of the die.
To obtain a desired extrudate shape, it is necessary to determine the corresponding die shape, which is an inverse problem. In industry, most of the present design procedures use a trial and error approach. A more efficient way to determine the corresponding die shape is to numerically predict the swell of the extrudate and solve the inverse problem to obtain the die corresponding to the desired extrudate shape. This inverse problem is studied extensively in recent years. Behr et al. [1] coupled finite element fluid flow simulations with a mathematical optimization algorithm to optimize the die shape. Gonçalves et al. [2] used an opti-Ellwood et al. [7] proposed a three-dimensional streamlined finite element method to predict the swelling of extrudates containing sharp edges. In this method, every element surface lies on a stream surface and every element edge lies on a streamline. The position of the free surface can then be obtained in two steps. The first step is to use that there is no mass penetration across the free surface and therefore the normal component of the velocity is zero. The second step is to use that a differential increment in the position vector defining the streamlines bounding a stream surface, is always parallel to the velocity vector along that streamline. Mitsoulis [8] validated this method for more complicated die shapes.
In the two methods discussed above, the positioning of the elements is very strict since all nodes have to be moved along streamlines or spines to obtain accurate free surface positions. Furthermore, these methods are steady state methods. This means that the swell of the extrudate can not be followed in time, but only the steady state dimensions of the extrudate can be obtained. We propose a new method, where we solve the temporal problem by only regarding the corner lines as material lines. This leaves the positioning of the elements free and allows for local mesh refinement. Using that the corner lines are material lines, equations have been derived from which the position of the lines in the two swell directions can be determined. Every line is solved as an additional problem to obtain the corner line positions. Every free surface separated by a corner line is treated as an individual problem on which a 2D height function is solved. Here the same discretization method introduced by Keunings is used [9] . The use of a height function to simulate extrudate swelling was already implemented in 2D by Choi and Hulsen [10] using an extended finite element method. The domain over which each separate height function is solved, is expanded using the obtained positions of the corner lines. This method is a transient technique, in contrast to the previous mentioned approaches. This means that it is possible to track the swelling in time making it a true transient method.
Wambersie and Crochet [11] developed a (pseudo)-transient 3D finite element method to predict extrudate swell for Newtonian fluids. An iterative calculation towards the steady state is performed with two different time steps: one for the velocity and pressure fields and one for the free surfaces, hence pseudo-transient. They treated the 2D free surfaces separated by a corner line as separate problems to obtain the free surface position. The difference with the method proposed here, is that they first calculated the positions of the free surfaces and than extrapolated the positions of the corner nodes, using the transverse movement of both surfaces, instead of describing the corner lines as material lines separately and using their positions to expand the free surface meshes in transverse direction. This means that in the here proposed method the positions of the corner lines are calculated instead of estimated and the 2D surface meshes are not fixed, but moving.
Fraggedakis et al. [12] proposed a new boundary fitted technique to describe free surface problems. The description of the mesh update is done by a simple set of elliptic equations. Like the ALE-method, this method accounts for smoothness, but also for orthogonality of the grid lines. This method was found to be accurate and robust, even for highly deformed boundaries.
The position where the fluid leaves the die is found to be highly singular. Comminal et al. [13] used a volume of fluid method together with a finite volume method to numerically calculate the planar extrudate swell for an Oldroyd-B model and a Giesekus model. They found that above a certain critical Weissenberg number the simulations were prone to surface instabilities, causing surface oscillations and wiggles on the free surface. Moreover they found that the calculations are very sensitive to mesh refinement.
The structure of this paper is as follows: first, we will present the proposed numerical method. The next step is to validate this model. We choose to do this with a 2D axisymmetric extrudate swell model. Therefore convergence of the 2D axisymmetric numerical method is tested first. Convergence results for a 2D axisymmetric problem will be shown for different constitutive models, to highlight several difficulties occur-ring due to the singular point at the die exit. The constitutive models that are used here are the Giesekus model, the linear Phan-Thien Tanner (PTT) model and the exponential PTT model. In order to better understand the difficulties in convergence, the rheological properties of the different models are compared. After the convergence results, validation of the 3D model will be done using the 2D axisymmetric results. Finally swell results will be presented for different complex die shapes for purely viscous fluids and for viscoelastic fluids. For the last case, the influence of different non-linear rheological parameters will be shown. Where possible, the results are compared to literature.

Problem description
The extrudate swell problem that is treated is that of a fluid, exiting a square die containing sharp edges. A schematic view of the problem is shown in Fig. 1 . The first part of the domain is the fluid contained in a die with half height H 0 . Here, a flow rate Q is applied at the inlet of the die. After a distance L die , the extrudate is modeled with a distance L from the end of the die.

Balance equations
It is assumed that the fluid is incompressible, inertia can be neglected and that there are no external body forces acting on the fluid. This leaves the following equations for the balance of momentum and balance of mass: where u is the fluid velocity and Ω is the fluid domain. Furthermore, is the Cauchy stress tensor, consisting of three contributions: where p is the pressure, I the unit tensor and 2 s D is the Newtonian (or viscous) stress tensor. Here, s is the Newtonian viscosity, = (∇ + ∇ )∕2 the rate of deformation tensor and represents the viscoelastic stress tensor.

Constitutive equations
The viscoelastic stress tensor for the models used here can be described as follows: where c is the conformation tensor and G the polymer modulus, which can be expressed as the ratio of the polymer viscosity p and the relaxation time , i.e. = p ∕ . The ratio between the solvent viscosity s and the zero shear viscosity 0 = s + p is called , i.e. = s ∕( s + p ) . The evolution of the conformation tensor c is given by: where D ()/ Dt denotes the material derivative and f ( c ) depends on the constitutive model used. In this paper three different constitutive models will be investigated: the Giesekus model, the linear Phan-Thien Tanner (PTT) model and the exponential PTT model. The Giesekus model is based on anisotropic relaxation, whereas the PTT model is a network model [14] . Moreover, the Giesekus model predicts a non-zero second normal stress difference N 2 .

Giesekus
The Giesekus model is given by: where is the mobility parameter that influences shear thinning.

PTT linear
The linear PTT model is given by: where is a non-linear parameter that is a measure for the amount of shear thinning.

PTT exponential
The exponential PTT model is given by: where is a non-linear parameter that is a measure for the amount of shear thinning. For the exponential model the relaxation time decays exponentially with stress for large stresses, meaning that at higher stresses the relaxation will be faster when compared to the Giesekus model and the linear PTT model.

Boundary and initial conditions
A schematic overview of the domain is shown in Fig. 1 . In order to be able to prescribe developed inflow conditions at the inlet of the boundary, first a subproblem of a channel with periodic boundaries is solved. A flow rate Q is prescribed to this channel as a constraint, using a Lagrange multiplier. The velocity and conformation tensor solution of this channel are prescribed as an essential boundary condition on the inlet boundary ( Γ in ) of the problem. At the walls of the die ( Γ die ) a no slip boundary condition is applied, whereas the tractions are prescribed to be zero at the free surfaces ( Γ free ). At the outlet ( Γ out ), there is only a velocity component in the x -direction, whereas the velocity components in the other directions are prescribed to be zero. The boundary conditions are given by: = 0 on Γ out , where u chan , c chan are obtained from the separate channel problem. The traction vector on the surface with outwardly directed normal n is denoted by = ⋅ . The free surfaces are material surfaces. Therefore, the following equation must hold: where ̇ is the velocity of the surfaces, u is the velocity of the material and n is the normal to the surface. Since a time derivative appears in Equation (5) for the viscoelastic fluid, an initial condition is needed for the conformation tensor c : which effectively means that the viscoelastic stress is initially assumed to be zero.

Arbitrary Lagrange-Eulerian formulation
The problem of extrudate swelling contains moving boundaries due to the movement of the free surfaces of the extrudate. We use a body fitted approach, which means that the mesh is not stationary. Therefore, the domain is described with a mesh that is moving in time in such a way that the mesh moves with the free surfaces, but not necessarily with the fluid. In order to do this, the governing equations are rewritten in Fig. 2. Movement of the corner lines (red) in an extrudate [4] . from a die containing sharp edges. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) the Arbitrary Lagrange-Eulerian (ALE) formulation [15] . Here, the grid is not stationary as in an Eulerian formulation but also does not follow material points as in a Lagrange formulation. The movement of the grid has to be taken into account in every equation that contains a material derivative, where it is split as follows: Here, ()∕ | | | denotes the time derivative in a fixed grid point and u m is the mesh velocity.

Description of free surfaces
In the case of a domain containing sharp edges, the corner lines (indicated in red in Fig. 2 ) will move in such a way that the sharp edge is maintained over a large distance in the extrudate. This means that the free surface can be treated as four separate free surfaces, connected by the corner lines, see Fig. 3 (a). Due to this connection the free surfaces will not only swell in y -direction, but also expand in z -direction due to the movement of the corner lines. Therefore, the projection on the coordinate planes of the swell surfaces is not constant, but expands according to the movement of the corner lines (indicated in blue in Fig. 3 (b)). The corner lines of the extrudate can be described as material lines (fluid particles remain on the line). This means that the motion in a plane normal to the tangential unit vector ̂ is material and therefore Equation (16) is enforced for any n normal to ̂ . The corner lines move in time in y -and z -direction and the y -and z -positions of every point on the line are also dependent on the x -position of that point, as is indicated in Fig. 4 . Therefore the position vector in every point of the line can be expressed as follows: From Fig. 4 it can be concluded that the velocity normal to the line in every point of the line can be obtained using the projection tensor: − ̂ ̂ . Equation (16) is obeyed by solving the following equation, using projection: where ̂ is the unit tangential vector, which can be expressed as follows: with 2 = 1 + ( ∕ ) 2 + ( ∕ ) 2 , ̇ is the velocity of the line and u is the velocity of the material and can be expressed as follows: Combining Eqs. (20) - (22) , gives the following set of equations (See Appendix A ): Therefore, the positions of the corner lines can be obtained by solving the following kinematic condition: where d is the position vector containing the positions in y -and zdirection = ( , ) , and u 2D is the velocity vector containing the velocities in y -and z -direction 2D = ( , ) . Therefore, four separate 1D line problems are created on which Equation (25) is solved for the problem in Fig. 2 .
In order to obtain the positions of the free surfaces a 2D height function is solved on the 2D free surface meshes [9] . If we take for example Γ free2 as shown in Fig. 1 , the height function is given by: where the subscript z indicates the direction of the expanding 2D ( x, z ) domain and u ml, z the corresponding mesh velocity. Four separate 2D surface problems are created for which Equation (26) is solved. For surfaces Γ free1 and Γ free4 , y and z are interchanged, due to the rotation of the surfaces with respect to Γ free2 . Once the positions of the corners lines are obtained, they can be used to expand the free surface meshes. The positions of the free surfaces in time can then be obtained by solving the 2D height functions compensated for the transverse movement of the corner lines perpendicular to the swell direction. This will cause the surfaces to expand and swell, as is shown in Fig. 5 .

Boundary and initial conditions for the free surface description
For the material line problems, the first point of the material line is attached to the die and does not swell: The initial condition on the height function d of the corner lines can be expressed as follows: with d 0 the position of the corner points of the die. For the free surface problems, the line of the free surface that is attached to the die does not swell: The initial condition on the height function h of the free surfaces can be expressed by: with H 0 the initial position of the free surface.

Numerical method
The finite element method is used to solve the governing equations. In order to solve the constitutive equation, the log-conformation representation [16] and SUPG are used for stability [17] . SUPG is also used for stability in the material line-and free surface problems. Furthermore, DEVSS-G is used for stability [18] .

Weak formulation
The weak formulation of the balance equations can be derived by multiplying the equations with test functions and integrating over the domain using partial integration and using Gauss' theorem. The weak form of this problem can now be formulated as follows: find u , p and c such that: ( for all admissable test functions v , q , H , . Furthermore, = (∇ + (∇ ) )∕2 , ( · , · ) denotes the inner product on domain Ω and and 1 are parameters due to DEVSS-G and SUPG stabilization respectively. Furthermore = log . More information on log-conformation stabilization and the function g can be found in [16] , whereas more information on the DEVSS-G method and the projected velocity gradient G can be found in [18] .

1D corner line problems
The weak form of the 1D corner line problems can now be formulated as follows: find d such that: for all admissable test functions . Here, 2 is a parameter for the SUPG stabilization.

2D free surface problems
The weak form of the 2D free surface problems can now be formulated as follows: find h such that: ( for all admissable test functions . Here, 3 is again a parameter for the SUPG stabilization.

Spatial discretization
For the velocity and the pressure isoparametric, tetrahedral P 2 P 1 (Taylor-Hood) elements are used, whereas for the conformation tetrahedral P 1 elements are used. For the 1D height functions of the corner lines quadratic line elements are used and for the 2D height functions of the free surfaces quadratic triangular elements are used. For the axisymmetric case, both quadratic and linear line elements are used for the 1D height function of the free surface. The mesh is generated using Gmsh [19] . Equations (34), (35) and (36) are solved using SUPG for stability. The SUPG parameters are obtained as follows: Here, is calculated in every integration point and h elem is the element size. For the conformation problem and the 2D free surface problems, h elem is defined using the method of Hughes et al. [20] . For the conformation problem SUPG = 1 . For the 1D corner line problems SUPG = 0 . 5 , h elem is defined by the length of the line element and | u | reduces to the velocity in x -direction u x . For the 2D free surface problems SUPG = 0 . 5 and | | = | surf | is the magnitude of the velocity in the plane of the free

Free surface problems
In Section 2.5 it was described that the domain of the 2D free surface meshes is expanded in z -direction. To do this, the ALE-method is used, where the mesh of the free surfaces moves with the material lines. In this way, the movement of the material lines is taken into account by moving the mesh nodes of the free surfaces with a mesh displacement d m . The mesh displacement can now be obtained by solving a Laplace equation [15] in order to guarantee a smooth mesh motion. From this mesh displacement, the mesh velocity can be obtained by numerically differentiating, using a second-order backward differencing scheme.

3D swell problem
The movement of the free surfaces is taken into account by moving the mesh nodes with a mesh displacement d m according to the movement of the free surfaces. This is done in the same way as for the free surface problems, by solving a Laplace equation. From this mesh displacement, the mesh velocity can be obtained by numerically differentiating, using a second-order backward differencing scheme.

Time integration
The numerical procedure for every time step will now be explained step by step. The positions of the free surfaces are obtained using a predictor-corrector scheme.
Step 1 : Update the position of the free surfaces, x free , in the bulk mesh using a prediction of the free surface positions. For the first time step the prediction of the positions equals the initial positions: free,pred = free, 0 . For subsequent time steps a second-order prediction of the free surface positions is used: Step 2 : Construct the ALE mesh by solving the Laplace equation. The new coordinates of the nodes are calculated with the obtained mesh displacement.
Step 3 : The mesh velocities can now be obtained by numerically differentiating the mesh displacement. In the first time step the mesh velocities are zero, since the height function is equal to the initial height H 0 . For subsequent time steps a second-order backward differencing scheme is used, using the updated mesh nodes: where Δt is the time step used.
Step 4 : Compute +1 and +1 . The method of D'Avino and Hulsen [21] for decoupling the momentum balance from the constitutive equation is applied. This means that at every time step the balance equations are solved using a prediction for the viscoelastic stress tensor, to find +1 and +1 .
Step 5 : After solving for the new velocities and pressures, the actual conformation tensor +1 is found using a second-order, semi-implicit extrapolated backward differencing scheme with conformation prediction for Equation (34) .  . 6. Schematic view of the 2D axisymmetric problem.
Step 6 : Update the position of the material lines by solving Equation (35) . For the first time step first-order time integration is used: whereas for subsequent time steps second-order time integration is used: Step 7 : Expand the surface meshes of the height function of the free surfaces in z -direction, using the displacement of the material lines. In this way it is known how the domain of the height function changes. Construct the ALE mesh by solving a Laplace equation. The new coordinates of the nodes are calculated with the obtained mesh displacement. The mesh velocities can now be obtained by numerically differentiating the mesh displacement. For the first time step a backward differencing scheme is used, using the updated mesh nodes. For, for example, Γ free2 this gives the following equation: where z m,surf is the displacement of the surface mesh in z -direction. For subsequent time steps a second-order backward differencing scheme is used: Step 8 : Update the position of the free surfaces by solving the evolution equation of the height function (36) . For the first time step first-order time integration is used, whereas for subsequent time steps second-order time integration is used, like explained in step 6.
Step 9 : Move the mesh nodes of the free surfaces of the bulk mesh, +1 free , to their time-accurate positions obtained by expanding the free surface meshes of the height functions in step 7 and solving the height functions in step 8.

Results
In this section the results of a convergence study of a 2D axisymmetric swell problem will be discussed to show what the difficulties in convergence are. This 2D axisymmetric problem will then be used to validate the new 3D corner-line method. Furthermore, swell results of different die shapes for the purely viscous case will be discussed as well as swell results for increasing Weissenberg number for the viscoelastic case.

2D axisymmetric problem
Convergence is studied for a 2D axisymmetric problem for the three different constitutive models. Later this 2D axisymmetric code will be used to validate the 3D method. For all models the same rheological parameters are chosen as shown in Table 1 . A schematic view of the 2D axisymmetric problem is shown in Fig. 6 .

Problem domain and meshes
Simulations are performed for an initial mesh with three elements along the radius and a 10 × refinement at the exit point at the upper wall of the die (See Fig. 7 ). This mesh is uniformly refined with a factor two up to six times for the convergence results. Table 2 shows the number of elements, number of elements along the radius at the die exit and the number of refinement steps with respect to the initial mesh M0 for the different meshes. In the simulations a Weissenberg number of 1 is used, calculated using the following definition: Fig. 7. Initial 2D axisymmetric mesh before uniform refinement steps.   Fig. 8 . The 1D height function for the 2D problem is given by: Applying the outflow boundary condition of Eqs. (12) - (14) , a shorter L / R ratio essentially means forcing the horizontal plateau in the steady state of the final extrudate shape to come earlier, since ℎ ∕ = 0 (steady state) and the boundary condition sets u y to zero in Equation (45) . This forces h / x to be zero. Fig. 8 shows that there is a slight difference in the final height of the extrudate at the outflow boundary for ∕ = 3 compared to longer extrudates. Since this difference is very small, we use ∕ = 3 to save computation time. This means, that essentially the swell of an extrudate that is solidified at an extrudate length of ∕ = 3 is modeled, which we believe is a relevant problem. All further simulation results use die ∕ = ∕ = 3 .

Die exit singularity
In Fig. 9 the pressure, scaled by the modulus G , is plotted for different element sizes over the upper vertical die wall and the free surface for the Giesekus model as indicated in Fig. 9 with the red line in the schematic view. It is clear that the exit point of the die is a singular point. Here a peak and a jump in the pressure occurs and it seems that this peak and jump in pressure increases for more refined meshes. When we zoom in on the area right after the die exit point we can clearly see that there remain oscillations in the pressure even for the most refined mesh (See Fig. 9 ). The jump in pressure was visible for all three models. The only observed difference was that the oscillations in the pressure are much smaller for the most refined mesh for the exponential PTT model as is shown in Fig. 10 . This figure also shows that the peak in pressure is less high for the exponential PTT model as for the Giesekus model.
Since the velocity has to change rapidly after the die, from zero to the uniform velocity at the free surface, the outer layer of the fluid is in extension [22] . Therefore it is interesting to compare the steady extensional viscosity of the three constitutive models. The result is shown in Fig. 11 . This figure shows that the exponential PTT model is the only model that predicts strain softening. This might explain why the peak and jump in the pressure is lower for the exponential PTT model and why the wiggles in the pressure seem to disappear for finer meshes. There are no big differences in rheological behaviour in shear for the three models. Fig. 12 (a) shows the volume of the extrudate for x ∈ [0, L ], divided by the initial volume over time, for the different meshes from Table 2 . This plot shows that the solution converges in the sense that upon mesh refinement the swell of all meshes tend to go to the same final volume change. Fig. 12 (b) is a zoomed-in version of Fig. 12 (a) and shows that the solution converges quite slowly. This was observed for all three constitutive models.

Mesh convergence
In order to investigate the convergence in a more quantitative manner, the following relative error is defined: where, V ( t ) is the volume of the extrudate at time t for a certain mesh,  Fig. 13 . This figure shows convergence of the final free surface position for the three different constitutive models upon mesh refinement. The convergence seems to be linear. This is lower than expected based on the order of interpolation of the elements. This is most likely caused by the singularity at the die exit, leading to a nonsmoothness of the solution. Furthermore, it was found that especially the Giesekus model and the linear PTT model are very sensitive to mesh refinement. When the mesh was not uniform refined, and therefore the elements at the die exit were too big, convergence was reached for the exponential PTT model but not for the Giesekus and linear PTT model.

Time convergence
In order to choose a suitable time step for the swell problem, time convergence is studied. Simulations are performed on mesh M3 for different time step sizes Δt for the three different constitutive models. For these simulations a Weissenberg number of Wi = 0 . 5 is used. There is a limit on how large the time step size can be due to the singularity at the die exit. This means that a smaller time step is required for more refined meshes, since the stress peak increases with decreasing element size, hence limiting the maximal time step size that can be used. This also means that for higher Weissenberg numbers smaller time steps are required. The time convergence is studied in a quantitative manner by calculating the error using Equation (46) for the volume of the extrudate at time = 5 calculated with different time step sizes Δt compared to a reference volume. Here, the reference volume is calculated on mesh M3 with a time step of Δ = 2 . 5 ⋅ 10 −3 . The result is shown in Fig. 14 . This figure shows linear convergence upon refining the time step. It also shows that the errors are quite small compared to the error in mesh convergence, indicating that it is only important that the time step is smaller than the time step limit due to the stress singularity at the die exit. Here, the convergence is again slower than expected based on the order of the time discretization scheme. This is most likely caused by the stress singularity at the die exit and the related non-smoothness of the solution.

Effect of solvent viscosity and model on surface instabilities
The value turned out to be important for stability. It was found that a relatively high value for is necessary to reach convergence for high Weissenberg numbers, i.e. Wi ≥ 1. For small Weissenberg numbers, instabilities occur when = 0 , especially for finer meshes. Wiggles in the free surface were observed and the simulations were terminated prematurely. The use of linear line elements to describe the height functions   Fig. 15 for the three models for linear and quadratic line elements. This figure shows that the exponential PTT model seems to be the most sensitive to , since wiggles appear for both linear and quadratic elements, whereas for the linear PTT model no wiggles appear for both, linear and quadratic elements. For the Giesekus model wiggles appear for quadratic elements, but not for linear elements. For the quadratic elements these wiggles seem to be convected further over the extrudate than for the exponential PTT model, before the code terminates. These instabilities arise from the stress singularity at the die exit and were also found for the Giesekus model by Comminal et al. [13] . The figures are zoomed in at the die exit, so not the whole extrudate is displayed, but only a portion of it. The wavelength of the wiggles is approximately 6 times the element size. For high Weissenberg numbers, the instability already occurs for a small value of ≈ 0.01 for the exponential PTT model. The stability results for different values of and Wi on mesh M4 are summarized in Table 3 . Similar surface instabilities as observed in our simulations are also observed during experiments and are called the sharkskin mechanism, where the surface defects arise due to high tensile stresses in the extrudate surface just after the die exit [23] . However, there is no evidence to prove that the surface instabilities occurring in the simulations represent this physical phenomenon.

Validation
The 3D corner-line method is validated using the results of the 2D axisymmetric model. For the sake of validation a 3D cylinder is mod- Fig. 11. Steady extensional viscosiy data for the three different constitutive models studied.

Table 3
Stability study for linear (L) and quadratic (Q) line elements on the free surface boundary for different values of and Wi on mesh M4. + indicates no problems, -indicates instabilities on the free surface and termination of the code and +/ -indicates that the simulation did not terminate prematurely, but instabilities appear on the free surface.
eled consisting of 4 free surfaces and therefore also of 4 material lines as is shown in Fig. 16 (a). The error of the final volume of the 3D cylinder with respect to the volume of the 2D axisymmetric reference mesh of Table 2 is calculated for different meshes using Equation (46) at time = 5 . Here, the 3D mesh is also locally refined near the die exit (See Fig. 16 (b)). The meshes used for the validation can be found in Table 4 . The validation is performed for the exponential PTT model. Fig. 17 shows the relative error of the 3D model and the 2D axisym-       metric model. The figure shows that both models converge to the same value upon mesh refinement, giving confidence that the 3D corner line method works.

3D viscous results
First the effect of the die shape for a purely viscous fluid is investigated. Where possible, the results are compared to literature. Simulations are performed on meshes with coarse element size 0.2 and refined element size 0.04. A time step Δ = 5 ⋅ 10 −3 * is used. Here, t * is a characteristic time scale defined as follows: * = 0 ∕ avg . In order to save computation time only a quarter of the domain is modeled, using appropriate symmetry boundary conditions. This is applicable to all die shapes used in this study. A schematic representation of the die shapes used in this study are shown Fig. 18 . Here, the black lines indicate the initial cross-section, the red dotted lines indicate the symmetry boundaries and the swell directions are indicated with letters.
Square: The method is now used to simulate the extrudate swelling of three different die shapes. In this section results for a square die are shown. Simulations are performed to calculate the swell percentages and the results are shown in Fig. 19 . The initial cross-section is indicated with the black line and the final simulated cross-section once the extrudate has reached steady state is indicated with the blue line. For the square die, the maximum swell indicated by the line OA in Fig. 18 (a) is 19%, and the minimum swell is the swell of the corners indicated by the line OB in Fig. 19    Rectangle: In practice, it is not common that extrusion dies are designed as a square. However, most flat extrusion dies containing sharp edges have a rectangular shape. Here, swelling of the extrudate from     a rectangular die with an aspect ratio of 2:1 is simulated. Simulations are performed to calculate the swell percentages. Results are shown in Fig. 21 . The maximum swell indicated by the line OC in Fig. 18 (b) is 25%, the swell indicated by the line OA is 12.35% and the swell indicated by the line OB is 4.6%. Fig. 21 shows that the material actually contracts at the corners. This was also concluded by Karagiannis et al. [6] . For the rectangle the material at the corner is pushed inward in order to maintain the force balance as the flow field changes.
The swell of a cross section of the rectangle at a distance ∕ 0 = 3 from the die in time is shown in Fig. 22 . Again, the swell in the corner stays behind compared to the maximum swell in the middle, and the sharp kink in the extrudate is maintained over a long time.
Cross: The third die shape that is investigated is that of a cross. The result of a swell simulation is shown in Fig. 23 . For these simulations the mesh is also refined near the sharp corner along swell-direction OB over the whole extrudate length. Here, the mesh size is the same as for the refinement near the die exit. Fig. 23 also shows a contraction of the corners in one direction as was also obtained from the rectangular die. The swell indicated by the line OA in Fig. 18 (c) is 15.1%, whereas the swell indicated by the line OB is 35.9% and the swell indicated by line OC is 7.5%. The swell of the corners points OC and OB are relatively high for the cross die compared to the other shapes. The swell of a cross section of the cross at a distance ∕ 0 = 3 from the die in time is shown in Fig. 24 . This Figure shows that first the corner points seem to be pushed upward compared to the swell of the upper horizontal line, due to the swelling of the sharp corner marked with line OB in Fig. 23 (a). All the extrudates also slightly shrink along the sharp edges when it first exits the die. This was also found by Trang-Cong and Phan-Thien

3D viscoelastic results
In this section the influence of the non-linear parameter of the Giesekus model and the exponential PTT model on the swell is investigated. Due to high computational costs, it was not possible to use more refined meshes than M2 in 3D. Since the error for this mesh was already very small (See Fig. 16 ) a mesh that is slightly coarser at the refinement point was chosen. Therefore, simulations are performed on a mesh with a coarse element size of 0.2 and a refined element size of 0.04 at the refinement point. A time step of Δ = 5 ⋅ 10 −3 is used. For high Weissenberg numbers and small non-linear parameter a smaller time step of Δ = 2 . 5 ⋅ 10 −3 is used. In all 3D viscoelastic simulations a value of = 0 . 09 is used. Giesekus: The steady state cross-sections of the extrudate from a square die for different Weissenberg numbers and different non-linear parameters are shown in Fig. 25 . This Figure shows that for higher Weissenberg numbers the swell increases. When the Weissenberg number is high enough, the sharp edges in the extrudate seem to disappear, as is shown for Wi = 2 and Wi = 3 in Fig. 25 (a). For increasing the swell decreases. Since increasing means more shear thinning, it can be concluded that shear thinning opposes swell.
PTT exponential: Fig. 26 shows the cross-sections of the extrudate from a square die for different Weissenberg numbers and different nonlinear parameters . Again it is shown that the swell increases for increasing Weissenberg number, and that the sharp edges eventually disappear if the Weissenberg number is high enough. Furthermore, the swell decreases with increasing non-linear parameter . Here, the reason is twofold. First of all for the PTT exponential model, the relaxation time is an exponential function of . This means that for high , is small and therefore the effect of elasticity gets less pronounced which leads to less extrudate swelling. Secondly, is a measure of shear thinning. A larger value for means that the fluid becomes more shear thinning, hence opposing the swelling mechanism.
Comparison: To compare the influence of the non-linear parameter and the Weissenberg number on the swell for the Giesekus and the exponential PTT model, the final change in volume of the extrudate in steady state is plotted against the Weissenberg number for different val-ues of the non-linear parameters. The result is shown in Fig. 27 . This figure shows that the Giesekus model predicts a higher swell value than the exponential PTT model. Furthermore, the swelling for small Weissenberg numbers seems to be smaller than the swelling of the purely viscous ( Wi = 0 ) case. This effect seems to get more pronounced when the non-linear parameter of the model (and therefore the amount of shear thinning in the material) is higher. It therefore seems that a small amount of elasticity can reduce the swelling compared to the purely viscous case. This was also found for the Giesekus model by Comminal et al. [13] in 2D.

Conclusions
In this paper a new transient 3D finite element method is developed to predict extrudate swell from complex die shapes containing sharp edges. In this method the sharp corner lines are described as separate material lines and the free surfaces as separate material surfaces. The positions of the corner lines can be obtained by solving a 1D height function. The obtained positions can than be used to expand the free surfaces. The positions of the free surfaces can than be obtained by solving a 2D height function.
First, convergence tests were performed with a 2D axisymmetric code to study the numerical difficulties of the problem and to be able to validate the new 3D method. Convergence tests were performed for the Giesekus model, the linear PTT model and the exponential PTT model and was found to be linear. This was lower than expected based on the order of interpolation of the elements and the order of the time discretization scheme. This is most likely caused by the stress singularity at the die exit and the related non-smoothness of the solution. Furthermore it was found that the Giesekus and linear PTT model were more sensitive to mesh refinement than the exponential PTT model. Wiggles remained in the pressure even for the finest mesh in the Giesekus and linear PTT model, while for the exponential PTT model these wiggles seemed to have disappeared. The peak and jump in the pressure kept on increasing with mesh refinement, but was found to be lower for the exponential PTT models than for the other two models. The value of the ratio of the polymer viscosity to the zero shear viscosity turned out to be important for stability. It was found that a relatively high value for was necessary to reach convergence for high Weissenberg numbers. The exponential PTT model seemed to be the most sensitive to the value of .
There is a clear difference in the rheological behaviour of the exponential PTT model compared to the other two models, since the exponential PTT model is the only model that predicts strain softening. This might explain the different behaviour of the exponential PTT model. It can therefore be concluded that the swell is highly dependent on the rheological properties of the fluid. Furthermore convergence remains challenging due to the stress singularity at the die exit.
Subsequently, the 2D axisymmetric code was used to validate the newly proposed 3D method. To this end, a cylinder with four swell surfaces, and therefore also four material lines was modeled and the results for different mesh sizes were compared to a 2D axisymmetric reference mesh. It was concluded that the final error in swell volume, compared to a 2D axisymmetric reference mesh converges to the same value for the 2D and 3D models.
After validation of the 3D method, simulations were performed for different complex die shapes for purely viscous fluids. The results compared favorably with results obtained from literature. Simulations for a square die for a viscoelastic fluid were performed for increasing Weissenberg number and increasing non-linear parameters of the Giesekus and exponential PTT model. It was found that the swell of the extrudate was greatly enhanced with increasing Weissenberg number and for high Weissenberg numbers the sharp edges eventually seemed to disappear. Increasing the non-linear parameter opposed the swell mechanism. Furthermore, the Giesekus model predicted a larger swell than the exponential PTT model and it was found that a small amount of elasticity seemed to reduce the swelling compared to the purely viscous case for both models. ∕ The third equation of (A.1) is given by: