# A level set method for simulating wrinkling of extruded viscoelastic sheets.

1 | INTRODUCTIONIn the last two decades, extrusion blow molding process became the leading technology to produce automotive fuel tanks of relatively large dimensions. [1] For such large parts, the extrusion is performed from annular dies or rectangular dies (eg, twin-sheet extrusion blow molding) of very large exit cross section aspect ratio. Experimental observations show that, in addition to the well-known die swell phenomenon, the extruded sheet (parison) has a wrinkled shape at its bottom edge portion composed of vertical bands.

Numerical modeling of wrinkle formation during polymer extrusion is almost inexistent in the scientific literature, most publications addressing the die swell simulation problem (see, eg, References [2-6]). Yousefi et al [2] and Yousefi and Atsbha [3] presented a comprehensive experimental study and numerical modeling of parison formation in extrusion blow molding. In their study, emphasis has been put specifically on die swell phenomenon due to stress relaxation and sag drawdown due to gravity. Tanoue et al [4] and Tanoue and Iemoto [5] carried out numerical simulation of parison formation process in blow molding. The flow field was divided into two regions, namely the extrudate swell region near the die lip and the parison formation region after the exit swell. Wrinkle formation has not been addressed. In the study related to the parison formation process by Tanifuji et al, [6] experimental observations showed that the extrudate parison has a wrinkled shape at its bottom portion formed in the initial stage of extrusion. Unfortunately, their numerical simulation could not predict such an irregular shape.

In order to elucidate the wrinkling phenomenon, Figure 1 shows a predicted shape image of a wrinkled high-density polyethylene (HDPE) sheet during extrusion. Experimental image (not shown for copyright protection) may be found in the work by Moitzheim. [1] Such wrinkles affect the overall shape of the extrudate, that is, diameter swell (for annular dies) or width swell (for slit dies) and may compromise subsequent steps of the manufacturing process. It is important for the design of the finished part to be able to control the shape and the dimensions of the extrudate along its entire length. Thus, the prediction of the occurrence of wrinkling during sheet extrusion is of great interest in the automotive industry.

It is well known that polymeric fluids may exhibit various extrudate distortions. [7] In particular, the melt fracture is observed when a polymer is extruded freely from a die at a rate exceeding a critical value. The diameter of the extrudate is no longer uniform and may exhibit various distortions, for example, sharkskin or helical shape. [7] With regard to the modeling of melt fracture instability, Atalik and Keunings'8' and Meulenbroek et al [9] conducted a nonlinear stability analysis for Poiseuille flow of viscoelastic fluids between plates and through a tube and confirmed a scenario in which Poiseuille flow of viscoelastic fluids exhibits a nonlinear subcritical instability due to normal stress effects. Meulenbroeck et al suggested that this nonlinear instability could be an important intrinsic route to melt fracture behavior in the absence of other mechanisms such as stick-slip phenomena. Grillet et al [10] also investigated numerically the linear stability of injection molding flow of an exponential Phan-Thien-Tanner (PTT) fluid to elucidate the mechanism of flow mark surface defects.

Wrinkling instability of freely extruded sheets studied in the present work is peculiar and differs from melt fracture instability and flow mark surface defects reported previously in the literature [9,10] as it is not a consequence of an unstable flow in the die. According to the work by Andeen, [11] extrusion of a sheet is an exception. The extruded materials desire to swell, while easily accommodated in the thickness direction, cannot be accommodated in the cross direction -without a wrinkling. In fact, wrinkling instability is a severe departure from the intended extrudate shape as illustrated in Figure 1 that could compromise subsequent steps of the manufacturing process.

A numerical study of wrinkling instability of freely extruded viscoelastic sheets appears to be absent from literature. Such an instability occurs particularly when the width-to-thickness ratio of the die exit section is large and the melt is viscoelastic, that is, the flow in the die is characterized by the existence of normal stress difference.

This study is intended to describe a 3D transient finite element method for predicting viscoelastic extrudate wrinkles that form during extrusion from rectangular dies involving large aspect ratio of the die exit cross section. In order to address complex free surface flow problems of a PTT viscoelastic fluid, the matrix logarithm-based formulation of the conformation tensor, [12] the stabilized discrete elastic viscous stress splitting method (DEVSS) [13] and the streamline upwind Petrov-Galerkin technique (SUPG) [14] are used. The single-phase level set method of Eulerian interface-capturing approach'15' is derived for viscoelastic fluid flow problems. Numerous factors that could influence the occurrence of wrinkling are investigated including the flow-rate, the material properties (Newtonian and viscoelastic fluids) and the die exit cross-sectional aspect ratio.

The rest of this article is organized as follows. First, the basic equations for viscoelastic fluid flows and the matrix logarithm of the conformation tensor in the context of PTT model are presented. Second, the single-phase level set formulation for free surface tracking is discussed. Third, Rayleigh's energy method is developed to determine the critical compressive stress for the onset of wrinkling in an extruded sheet. Finally, numerical analyses of free sheet extrusion are conducted to determine the critical conditions for onset of wrinkling.

2 | BASIC EQUATIONS FOR VISCOELASTIC FLUID FLOWS

We consider the time-dependent isothermal flow of an incompressible viscoelastic fluid. Let us denote the domain containing the viscoelastic liquid by [[OMEGA].sub.L] (gray region in Figure 2), and the domain boundary between the viscoelastic liquid and the surrounding air by [[GAMMA].sub.L]. Expressed in Eulerian form, the general conservation laws for linear momentum and mass in the liquid region read

-[nabla]p + [nabla] - [tau] + [rho]g = p [[partial derivative]u/[partial derivative]t + (u - [nabla]), (1)

[tau] = [[tau].sub.s] + [[tau].sub.p] = 2[[eta].sub.s]D + [[tau].sub.p], (2)

[nabla] * u = 0 (3)

Here, u is the fluid velocity, [rho] is the fluid density, g is the gravity and p is the hydrodynamic pressure. In general, the extra stress tensor, [tau], can be split into a polymeric contribution [[tau].sub.p] and a solvent contribution [[tau].sub.s] with [[eta].sub.s] the solvent viscosity and D the rate-of-deformation tensor.

We emphasize that the presence of a purely viscous component has a great impact on the mathematical nature of the full set of governing equations and also extends the domain of convergence of the mathematical method in a significant way. [16,17] Keunings [17] showed that viscoelastic fluid models without a purely viscous component can exhibit a variety of hyperbolic phenomena, including change of type of the governing equations and propagation of wave. As in the present study, we are mainly interested in wrinkling instability of polymer melts (without solvent), the purely viscous component will not be included ([[tau].sub.s] = 0) to avoid any alteration of the solution. Therefore, Equation (2) reduces to [tau] = [[tau].sub.p]. In the sequel, only the tensor [tau] will be used.

Boundary conditions (velocity components or surface force components along the boundary of the flow domain) have to be specified:

u = [bar.u], at the die inlet boundary, (4)

u = 0, on the die wall. (5)

We further assume that no external forces act on the free surface [[GAMMA].sub.L], that is,

t = [tau] * n - pn = 0. (6)

In the present study, we wish to solve a 3D viscoelastic flow in which normal stress effects are of interest. For an isothermal flow, the simplified PTT rheological constitutive model is given by the following equation [18]:

A([tau])[tau] + [lambda][??] =2[[eta].sub.p]D, (7)

where 2 is a relaxation time, [[eta].sub.p] is a viscosity coefficient, A is a model-dependent scalar function, and [??] denotes Oldroyd1S upper convected derivative. For a PTT fluid, the scalar function A can take the form of the exponential equation,

A([tau]) = exp ([epsilon][lambda]/([[eta].sub.p] tr([tau]), (8)

where [epsilon] is a nondimensional scalar parameter.

The constitutive equation, Equation (7), is completed by specifying the relationship between the polymeric stress contribution [tau] and the conformation tensor c. This relationship is taken to be of the form

[tau] = [[eta].sub.p]/[lambda](c-Z). (9)

In terms of the conformation tensor, c, the PTT constitutive model, Equation (7), can be written as

[??] = [nabla] [u.sup.T] * c + c * [nabla]u + f(c), (10)

where the tensor function, f, is given by

f(c)= - 1/[lambda] |exp[[epsilon](tr(c)-3)](c - I). (11)

In order to complete the mathematical description of the viscoelastic fluid flow problem, we must specify appropriate boundary conditions. The nature of those boundary conditions is intimately connected to the mathematical nature of the governing equation. As Equation (10) is hyperbolic, it requires an initial distribution of the viscoelastic stress tensor in the filled domain together with boundary conditions on a line or boundary crossing the characteristics of the first-order system of Equation (10), that is, at the die inlet section. In the case of viscoelastic fluid flow problems, the usual assumption of a fully developed flow at the inlet section allows the computation of the corresponding extra-stress components, which are then applied as Dirichlet conditions.

3 | LOGARITHM OF THE CONFORMATION TENSOR

We follow the method developed in References [12,19] to derive the evolution equation of the conformation tensor c in logarithm base. In that method, the diagonalizing transformation of c can be written as

c = R[??][R.sup.T] (12)

where [??] is a diagonal matrix whose diagonal elements, [[??].sub.kk], are eigenvalues of c and R is an orthogonal tensor whose columns are eigenvectors of c. Therefore, the matrix logarithm of the conformation tensor can be written as

[psi] = logc = R(log[??])[R.sup.T]. (13)

After some mathematical manipulations, the transformation from Equation (10) to an evolution equation for [psi] leads to

[[??].sub.ij] = [R.sub.ik][[??].sub.km][R.sub.jm], (14)

in which, the components of the tensor [bar.N] are defined as

[mathematical expression not reproducible] (15)

and

[mathematical expression not reproducible] (16)

Here, L = [nabla] [u.sup.T] is the transpose of velocity gradient tensor. When [[??].sub.mm] [right arrow] [[??].sub.kk], we get from Equation (16):

[mathematical expression not reproducible] (17)

4 | FREE SURFACE TRACKING: SINGLE-PHASE LEVEL SET FORMULATION

Unlike the conventional level set method, with the single-phase level set-based method, the momentum, continuity, and viscoelastic equations are only computed in the liquid phase of the fluid domain and the level set function is used as a tracking device to locate the actual position of the free surface on a fixed grid. [15]

For the sake of clarity, let us first recall the conventional level set method. Assume that [[OMEGA].sub.L] is the liquid domain, [[OMEGA].sub.a] is the air domain of the full domain [OMEGA] (Figure 2). As the liquid domain evolves over time, the following notation [[OMEGA].sub.L] should be used to define the fluid domain at time [t.sub.n]. to simplify the notation, [[OMEGA].sub.L] is still used instead of [[OMEGA].sup.n.sub.L]. The interface [[GAMMA].sub.L] between the two subdomains is defined by the level set scalar function [phi] that takes the form of a signed distance to the interface [[GAMMA].sub.L] and satisfies the following conditions [20]:

[mathematical expression not reproducible] (18)

Since the interface, [[GAMMA].sub.L], moves with the liquid, its motion results from the following advection equation in the Eulerian coordinates

[partial derivative][phi]/ [partial derivative]t + u * [nabla][phi] = 0, [phi](x,t = O) = [[phi].sub.0](x). (19)

The level set function [phi] is conveyed using the velocity field u provided by the solution of the momentum and the continuity equations, Equations (1) to (3). This function can quickly cease to be a signed distance function. [21,23] To remedy the problem, a re-initialization technique as proposed by Sussman et al [22] has to be applied to ensure that the scheme accurately conserves mass. The re-initialization benefit from the fact that a distance function must satisfy [mathematical expression not reproducible], where [mathematical expression not reproducible] is the Euclidean norm.

In this study, however, we follow the method proposed by Ville et al [24] and Bonito et al. [25] In that method, a hyperbolic tangent filter is applied to the distance function to get a smooth truncation far from the interface, that is,

[??](x,t) = [beta]tanh (d(x,t)/[beta]), (20)

where [beta] > 0 is the width of the hyperbolic tangent filter and will be of the order of the mesh size in our numerical simulation. It will be automatically calculated.

Also, the re-initialization is embedded in the transport equation model, avoiding it to be a separate step during calculation. As in our computation, we implemented an "on the fly" re-initialization algorithm proposed in References [24,25], which consists of replacing Equation (19) by the convected re-initialization equation as (to simplify the notation, [phi] is still used instead of [??])

[mathematical expression not reproducible] (21)

where the signed function sign ([phi]) is defined by [25]

[mathematical expression not reproducible] (22)

with [C.sub.[phi]] > 0.

In Equation (21), the parameter [theta] > 0, called the re-initialization relative speed, is defined as [24,25]

[mathematical expression not reproducible] (23)

where [mathematical expression not reproducible] is the infinity norm and [C.sub.[theta]] > 0.

In the single-phase level set method, the computational domain used for solving the governing equations, Equations (1) to (3), is defined to be the union of all liquid elements. An element (tetrahedron) of the mesh is said to be liquid if, at least, one of its vertexes has a distance function [phi] > 0. [15,26] We note that outside the fluid region, u = 0 and Equation (21) reduces to

[mathematical expression not reproducible] (24)

Due to the hyperbolic nature of Equation (24), the solution evolves as a distance function close to the zero level set first, then it propagates outward in the whole domain as

[phi](x,t) = [beta]tanh(d(x,t)/[beta]) (25)

Note that the solution of Equation (24) is such that [mathematical expression not reproducible] in any small neighborhood of [phi] = 0. [25]

5 | SPACE-TIME DISCRETIZATION AND SOLUTION PROCEDURE

We follow the method developed by Kabanemi and Hetu [27] for solving viscoelastic fluid flow problems. At each time step, the resolution of the full model is performed in a decoupled fashion using an iterative scheme. The flow domain is discretized by means of tetrahedral finite elements. Implicit Euler approximation to the time derivative in the momentum and constitutive equations is used. In order to solve the momentum and the continuity equations in the fluid domain [[OMEGA].sub.L], the stabilized discrete elastic viscous stress splitting (DEVSS) method, developed by Guenette and Fortin1131 is used. In that method, an additional unknown, the rate of deformation tensor D = ([nabla][u.sup.T] + [nabla] u)/2 is introduced for stability purposes and computed from u in a straightforward manner using a finite element approximation of D in the least-squares sense, that is,

[mathematical expression not reproducible] (26)

where the superscript T denotes the transpose operation and [S.sub.D] denotes the function space defined on [[OMEGA].sub.L] for D. Having calculated D, the new velocity [u.sup.n] and pressure [p.sup.n] at current time step, [DELTA][t.sub.n] = [t.sub.n] - [t.sub.n-1], can be computed by solving the Stokes problem, with the viscoelastic stress [tau] treated as a given body force. In this study, according to the general theory developed by Hughes et al [28] and Behr and Tezduyar [29] we use a formulation with the stabilization terms from the square residuals of the continuity and momentum equations. It is stated as follows: find ([u.sup.n], [p.sup.n]) [member of] [S.sub.u] X [S.sub.p] such that

[mathematical expression not reproducible] (27)

[mathematical expression not reproducible] (28)

where [[eta].sub.a] is an artificial viscosity to achieve numerical stability, Su and Sp denote the function spaces defined on [[OMEGA].sub.L] for u and p, respectively. The stabilization parameters [[alpha].sub.c] and [[alpha].sub.m] follow standard definition. [29] Such a stabilized method allows the use of velocity and pressure interpolants that do not satisfy the Babuska-Brezzi condition such as the linear equal order interpolation functions. The classical Courant-Friedrichs-Lewy (CFL) criterion is used for the automatic selection of the time step [DELTA][t.sub.n]. Its value is limited by the smallest element size in the fluid domain and the highest velocity. Unconditional stability is guaranteed with this method.130'

Once the velocity u and the rate-of-deformation tensor D are computed, the constitutive equation is integrated in the fluid domain [[OMEGA].sub.L] with known kinematics by using the streamline upwind Petrov-Galerkin technique (SUPG) as developed by Carew et al, [14] for the matrix logarithm tensor [psi]. The conformation tensor c is then calculated from

c = exp([psi]) = R[exp([??])][R.sup.T]. (29)

A fixed-point algorithm is used to iterate between the solution of the momentum and constitutive equations at each time step, thereby making the overall algorithm coupled in a segregated manner. The resulting formulation is well posed.

In our numerical simulation, the relative convergence criterion is [10.sup.-4], that is,

[mathematical expression not reproducible] (30)

where [u.sup.p] and [[tau].sub.p] are the velocity and stress from the previous Picard iteration, respectively. Generally, the iterations are interrupted when the error for all variables lies below [10.sup.-3] or [10.sup.-4]. [31]

The level set equation is solved in the whole computational domain. Finally, a newly filled domain is obtained for the next time step. The methodology used in the present work is summarized in Reference [27],

The numerical methods just described have been implemented into dFEMwork, NRCs finite element software toolkit for large scale parallel computing. This toolkit provides distributed data structures such as matrices and vectors, parallel preconditioned Krylov iterative solvers and dynamic data redistribution. Domain decomposition was performed using ParMETIS parallel graph partitioner. [32] Algebraic systems were solved using incomplete factorization ILU(O) conjugate gradient iterative method.

6 | RAYLEIGH'S ENERGY METHOD FOR WRINKLING ANALYSIS

In order to determine the critical condition for wrinkling of extruded sheets in an approximate way, we use the energy method (see, eg, Reference [33]). Rayleigh's method is based on the principle of conservation of energy. The potential energy in our system includes strain energy that is proportional to elastic deformations and the work done by the applied forces.

Consider an elastic extruded sheet with a rectangular cross section of thickness h and width L as shown in Figure 3. The sheet is flat (that is, in x-z plane, Figure 3A), the extrusion is along the z-axis while x-axis is along the width.

In applying the Rayleigh' energy method to the wrinkling of a sheet, for which v(x) is the transverse displacement (that is, the displacement in y direction or the out-ofplane displacement along the mid-section of the sheet) at section x (see Figure 3B), one writes the strain energy per unit length along z, stored in the wrinkling (bending) sheet as

[mathematical expression not reproducible] (31)

where M = D/R is the bending moment per unit length along z, R is the radius of curvature of the deformed sheet in y direction, 1/R = [[partial derivative].sup.2]v/[partial derivative][x.sup.2] is the curvature, D is the flexural rigidity of the sheet and is given by

D = [Eh.sup.3]/ 12(1 - [v.sup.2]), (32)

where E is the Young's modulus of the sheet and v is the Poisson's ratio.

Let [DELTA]l be a virtual displacement along x-axis compatible with boundary conditions. The work of external forces F corresponding to the virtual displacement, that is, the energy causing wrinkle formation or bending of the sheet is

[mathematical expression not reproducible] (33)

where F is the applied load per unit length along z-axis. The potential energy, U, associated with virtual displacement is the difference

U = W - V. (34)

According to the principle of minimum potential energy, the system is in a stable equilibrium state if, during a virtual displacement superimposed on the equilibrium state, U is positive. The equilibrium is unstable if U [less than or equal to] 0. The critical wrinkling load is obtained when

V = W (35)

The critical load depends on the choice of the transverse displacement solution. Assume a sinusoidal wrinkle solution function, with amplitude A and wavelength 2 L of the form

v(x) =Asin [pi]x/L (36)

Thus, the critical load per unit length along z-axis is

[F.sub.c] = [Eh.sup.3]/12(1 - [v.sup.2]) [([pi]/L).sup.2], (37)

Equivalently, the critical stress to cause wrinkling is

[[sigma].sub.c] = [F.sub.c]/h = [[pi].sup.2]E/12(1 - [v.sup.2]) 1/[chi square], (38)

where [chi] = L/h is the width-to-thickness aspect ratio of the die exit cross section. It follows from Equation (38), which represents an approximate model for wrinkling analysis, that the critical compressive stress, [[sigma].sub.c], for the onset of wrinkling of an elastic sheet scales like [[sigma].sub.c]~1/[chi square], with a significant drop for [chi] [much greater than] 1. Therefore, for very large [chi], as soon as a slight compressive stress in the width direction is induced in the extruded sheet next to the die exit lip, it could potentially cause wrinkling.

7 | NUMERICAL RESULTS

The practical utility and effectiveness of the proposed numerical scheme is demonstrated by solving a fully three-dimensional die swell problem. The flow is characterized by the Reynolds number defined as, Re = p[bar.U]h/[eta], ranging from 1 * [10.sup.-7] to 7 * [10.sup.-7], where p is the density, U is the average velocity of the flow in the die, h is the thickness of the die, and [eta] is the viscosity. Although Re is small, the inertia terms cannot be neglected in our analysis. Values of the numerical parameters for the level set equations are [C.sub.[phi]] = [10.sup.-4] (Equation (22)) and [C.sub.[theta]] = [10.sup.-2] (Equation (23)). Simulations were run on a Cluster of UNIX PCs, using 64, 2.67-GHz-processors.

7.1 | Extrusion from a rectangular cross-sectional die of large aspect ratio: Newtonian case

In order to validate the numerical method developed in this study, we consider the extrusion from a rectangular die of large width-to-thickness ratio, that is, [chi] = L/h = 62.5 (thickness h = 8 mm and width L = 500 mm). The sketch of the problem is illustrated in Figure 2. The flow is assumed to be creeping. The main flow is along the z-axis. We impose noslip conditions on the walls of the die and vanishing surface forces on the free surface. We assume that the die is long enough, that is, [L.sub.d] = 125 h (die thickness) = 1000 mm, so that we may impose a bi-quadratic velocity profile at the die inlet cross section, which evolves to a fully developed profile before emerging from the die. The inlet velocity profile has the following form:

w(x,y) = [V.sub.max] [[(L/2).sup.2] - [x.sup.2]][[(h/2).sup.2] - [y.sup.2]]/ [(L/2).sup.2][(h/2).sup.2] (39)

where [V.sub.max] = 15 mm/s. We wish to predict the thickness swelling ratio defined by

[S.sub.h] = [h.sub.e]/h (40)

where h is the thickness of the die exit cross section and [h.sub.e] is the thickness of the extrudate cross section. We also define the width swelling ratio as

[S.sub.w] = [L.sub.e]/L, (41)

where L is the width of the die exit cross section and [L.sub.e] is the width of the extrudate cross section.

In this simulation, the value of the Newtonian viscosity is that of the zero shear rate viscosity of HDPE given in Table 1 and its density is [rho] = 746 kg/[m.sup.3]. The artificial viscosity [[eta].sub.a] in the DEVSS method is the same as that of the zero shear rate viscosity of HDPE.

Finally, we emphasize that we have solved the problem on a relatively fine finite element mesh (not shown here) consisting of 15.7 * [10.sup.6] tetrahedrons and 2.7 * [10.sup.6] nodes. The effectiveness of this mesh is demonstrated in the following results.

The computed die swell results are reported in Figure 4. The thickness swelling ratio (Figure 4B) is [S.sub.h] = 1.19, which only slightly differs from the theoretical swelling ratio of 1.2 obtained by Tanner [34] Such a result gives strong evidence that the mesh density used is large enough to provide accurate and converged solution with respect to the mesh refinement. The width swelling ratio at a distance of z = 100 mm from the die exit cross section is [S.sub.w] = 1.032 (Figure 4D). In Figure 5, we show the effect of sag on the shape of the extruded sheet after 80 seconds of extrusion. The extruded sheet length was increased by a factor of 8% due to the sag effect.

7.2 | Extrusion through a rectangular cross-sectional die of large aspect ratio: Viscoelastic case

We now consider the extrusion of a viscoelastic PTT fluid from a rectangular die of large width-to-thickness ratio. The sketch of the problem is shown in Figure 2. The polymer investigated is the HDPE Lupolen 4261A. The values of the parameters used in the simulation are summarized in Table 1. The parameter [epsilon] of the PTT model has been determined as best fitting of the shear viscosity data. In Figure 6, we have plotted the shear viscosity of the PTT model together with the HDPE's experimental measurements. The agreement is reasonable within engineering approximations. The HDPE Lupolen 4261A exhibits high shear-thinning behavior in the shear rate range of the extrusion process. Also shown in the same figure is the first normal stress difference [N.sub.1]/2.

We assume that the die is long enough, that is, [L.sub.d] = 125 h (die thickness) = 1000 mm, so that we may impose a bi-quadratic velocity and extra-stress profiles at the die inlet cross section, which evolves to a fully developed profile before emerging from the die. The inlet velocity profile is given by Equation (39). We also impose no-slip conditions on the walls of the die and vanishing surface forces on the free surface. In order to characterize the elasticity of the polymer, we define the Weissenberg number Wi by

Wi = [lambda][[??].sub.w], (42)

where [lambda] is the relaxation time, and [[??].sub.w] is the shear rate at wall in the upstream fully developed flow. [35]

Let us first consider the extrusion from a rectangular die exit cross section with an aspect ratio [chi] of 12.5 (that is, thickness h = 8 mm and width L = 100 mm). The inlet velocity profile is given by Equation (39), with [V.sub.max] = 5 mm/s, that corresponds to a Weissenberg number of about Wi = 40. The gravity is also included in the simulation. We show in Figure 7 the predicted parison shape at various times during extrusion. Throughout the extrusion, the sag due to the parison weight becomes prominent (time = 333 seconds), and a swelling reduction is observed at mid-length of the extrudate. In Figure 8, we show the extrudate swell in the cross section, at z = 60 mm from the die exit. The thickness swelling ratio is 1.45, while the width swelling ratio is 1.34. It is worth to mention that these swelling ratio are not uniform along the extrudate due to the sag. The distribution of the first normal stress difference [N.sub.1] = [[tau].sub.zz] - [[tau].sub.yy] along the central line, in the flow axis, is presented in Figure 9A (the flow is along the z-axis while the shear is oriented with the y-axis). It starts from zero in the die, increases to reach a maximum just before the exit and then decreases to become negative just before the die exit. A minimum negative value of [N.sub.1] - 0.065 MPa is found just after the die exit. By applying the energy method for wrinkling (see Equation (38)), the critical compressive stress for the sheet to wrinkle is [[sigma].sub.c] = -1.317 MPa (the aspect ratio of the die is [chi] = 12.5, the Young's modulus E for the HDPE is 220 MPaJ361 and the Poisson's ratio is v = 0.35). As the minimum compressive stress induced in the extruded sheet soon after exit is high (-0.065 MPa) compared to the critical compressive stress [[sigma].sub.c] for the sheet to wrinkle, that is, [N.sub.1] > [[sigma].sub.c] (note that we are comparing negative stresses), the extrudate does not wrinkle, in accordance with the energy method presented previously.

To understand the effect of the aspect ratio [chi] of the die exit cross section, we consider the case where its value was increased from [chi] = 12.5 to [chi] = 62.5 by multiplying the width by a factor of 5 (that is, thickness h = 8 mm and width L = 500 mm). In this case, the energy method developed previously (see Equation (38)) shows that the critical compressive stress for wrinkling increases drastically from [[sigma].sub.c] = - 1.317 MPa for [chi] = 12.5 to [[sigma].sub.c] = - 0.0524 MPa for [chi] = 62.5. In fact, the critical compressive stress for wrinkling scales like 1//2. In other words, an extruded sheet of very large aspect ratio provides lower critical stress and constraints for wrinkling.

Let us study the shape of the extruded sheet at Weissenberg number, Wi, below critical value, for a fixed aspect ratio [chi] = 62.5. In Figure 10, we show the extrudate shape after 160 seconds of extrusion at Wi = 27.8 (corresponding to [V.sub.max] = 3.5 mm/s in Equation (39)), just below the critical Weissenberg number for the onset of wrinkling. In the same figure, we also show the cross-sectional shape of the extruded sheet at a distance of 60 mm from the die exit. The distribution of the first normal stress difference [N.sub.1] = [[tau].sub.zz] - [[tau].sub.yy] along the central line axis, in the flow direction, is presented in Figure 9B. A minimum negative value of [N.sub.1] = -0.0512 MPa is found just after the die exit. As this compressive stress induced in the extruded sheet soon after exit is high ([N.sub.1] = -0.0512 MPa) compared to the critical compressive stress ([[sigma].sub.c] = -0.0524 MPa) for the sheet to wrinkle, that is, [N.sub.1] > [[sigma].sub.c], the extrudate remains flat (see Figure 10), in accordance with the energy method presented previously.

Let us now study the effect of increasing the Weissenberg number close to the threshold of the instability, on the shape of the extruded sheet for a fixed aspect ratio [chi] = 62.5. In Figure 11, we show the extrudate shape at Wi = 29 (corresponding to [V.sub.max]= 4 mm/s in Equation (39)). Unlike the previous results at Wi = 27.8, the extruded sheet no longer remains flat. The crosssectional shape of the extruded sheet shown in Figure 11 deforms slightly out of the plane half its width. We conclude that the value of Wi = 29 corresponds to the critical Weissenberg number for the onset of wrinkling for an extruded sheet with aspect ratio [chi] = 62.5. The distribution of the first normal stress difference [N.sub.1] = [[tau].sub.zz] - [[tau].sub.yy] along the central line axis, in the flow direction, is presented in Figure 9C. A minimum negative value of [N.sub.1] = -0.0565 MPa is found just after the die exit, while the critical compressive stress for the sheet to wrinkle is [[sigma].sub.c] = -0.0524 MPa (from energy method). As the compressive stress induced in the sheet next to the die exit is slightly below the critical compressive stress for the extruded sheet to wrinkle, that is, [N.sub.1] < [[sigma].sub.c], wrinkles of very low amplitude occur, as highlighted in Figure 11.

By increasing further the Weissenberg number Wi to 40 (corresponding to [V.sub.max] = 5.5 mm/s in Equation (39)), sinusoidal wrinkles of high amplitude develop soon after the sheet is extruded from the die exit as shown in Figure 12.

The distribution of the first normal stress difference [N.sub.1] along the central line, in the flow axis, is presented in Figure 9D. It starts from zero in the die, increases to reach a maximum just before the die exit, then decreases to become negative and reach a minimum negative value just after the die exit. After a certain distance from the die exit, the normal and shear stresses relax to zero in absence of confinements. Therefore, soon after extrusion from the die exit, the extruded sheet behaves as a rectangular sheet subjected to a compressive residual stress in the width direction. Upon wrinkling, the sheet deforms to relax compressive stresses. The magnitude of such stresses that cause wrinkling depends solely on the magnitude of the normal stress difference developed in the die. The minimum compressive stress induced in the extruded sheet soon after exit from the die is [N.sub.1] = -0.065 MPa (see Figure 9D). As this compressive stress is below the critical compressive stress ([[sigma].sub.c] = -0.0524 MPa) for the extruded sheet to wrinkle, that is, [N.sub.1] < [[sigma].sub.c] (note that we are comparing negative stresses), wrinkling occurs, in accordance with the energy method presented previously.

Predicted extrudate wrinkles at various times during extrusion are shown in detail in Figure 13 (front view) and in Figure 14 (perspective view). By comparing these results with the previous ones (see Figure 7) obtained at same Weissenberg number (Wi = 40) but with a smaller die exit cross-sectional aspect ratio ([chi] = 12.5), we observe the absence of wrinkles for a smaller aspect ratio. All these findings are consistent with the energy method which shows that the critical compressive stress [[sigma].sub.c] for sheet wrinkling scales like l/[chi square].

In the particular case of extrusion of a Newtonian fluid from a rectangular die with a large aspect ratio [chi] = 62.5, shown in Figures 4 and 5, the absence of normal stress difference in the die had the effect of maintaining a flat extruded sheet without wrinkling. In fact, in that case, there is no driving force to trigger the formation of wrinkles in the extruded sheet. We therefore conclude that wrinkling instability is purely elastic.

Finally, various wrinkle patterns may be predicted, depending on the aspect ratio [chi] and the flow rate. In Figure 15, we summarize the results for [chi] = 62.5 and two inlet velocities of [V.sub.max] = 5.5 mm/s in Equation (39) (corresponding to Wi = 40) and [V.sub.max] = 7 mm/s (corresponding to Wi = 52). Wrinkle patterns in the form of sine-like and cosine-like are, respectively, predicted. The amplitude of wrinkles increases with increasing flow rate.

8 | CONCLUSIONS

The objective of the current research was to investigate, by means of transient finite element method, wrinkling of extruded viscoelastic sheets with large width-to-thickness ratio [chi] of the die exit cross section. In order to be able to predict the critical conditions for an extruded viscoelastic sheet to wrinkle, we combined the matrix-logarithm-based formulation of the conformation tensor for a PTT viscoelastic model and the single-phase level set method based on a hyperbolic tangent filter. The effects of die exit cross-sectional aspect ratio and flow rate or Weissenberg number on wrinkling of extruded sheet were investigated in detail. For a given sheet cross-sectional aspect ratio and material properties, an approximate model based on Rayleigh's energy method shows that the critical compressive elastic stress, [[sigma].sub.c], for the onset of wrinkling of an elastic sheet scales like [[sigma].sub.c]~1/[chi square]. Numerical results based on transient finite element method correctly predicted the occurrence of wrinkle formation during sheet extrusion process in accordance with the threshold value of the critical compressive stress for wrinkling predicted by the energy method. Wrinkling instability studied in this work differs from melt fracture instability and flow mark surface defects reported previously in the literature as it is not a consequence of an unstable flow in the die, but rather induced by the combined effects of normal stress difference developed in the die and large die exit cross-sectional aspect ratio, which are prerequisites for wrinkling. In the particular case of extrusion of a Newtonian fluid from a rectangular die with a large die exit cross-sectional aspect ratio, the absence of normal stress difference in the die had the effect of maintaining a flat extruded sheet without wrinkling. In that case, there is no driving force to trigger the formation of wrinkles in the extruded sheet. It follows that the elastic properties of the polymer play a major role in the occurrence of wrinkle instability in extruded sheets. Such findings could be subsequently used in the die design optimization.

DOI: 10.1002/pen.25409

ACKNOWLEDGMENTS

The authors would like to acknowledge the members of NRCs SIGBLOW Consortium and the Advanced Manufacturing Program (AMP) for their support and input in our R&D activities. The authors would like also to thank the members of the Simulation and Numerical Modeling team for their valuable assistance in this project.

AUTHOR CONTRIBUTIONS

Kalonji K. Kabanemi is a project leader who developed and implemented viscoelastic models and algorithms for 3D die swell and wrinkling problems. He also wrote the present draft. Jean-Philippe Marcotte participated in the implementation of models and algorithms. He also contributed in running simulations, collecting numerical results, and writing the present draft.

ORCID

Kalonji K. Kabanemi [ID] https://orcid.org/0000-0003-3995391X

Received: 3 February 2020 | Revised: 4 April 2020 | Accepted: 22 April 2020

Kalonji K. Kabanemi [ID] | Jean-Philippe Marcotte

National Research Council of Canada, Boucherville, Quebec, Canada

Correspondence

Kalonji K. Kabanemi, National Research Council of Canada 75 De Mortagne Blvd. Boucherville, Quebec, J4B 6Y4, Canada. Email: kalonji.kabanemi@nrc.ca

Reproduced with the permission of the National Research Council Canada.

REFERENCES

[1] J. Moitzheim, Fuel Systems for LEVIII and PZEV Vehicles, Kautex Maschinenbau, Germany 2012, p. 15.

[2] A. M. Yousefi, P. Collins, S. Chang, R. W. DiRaddo, Polym. Eng. Set. 2007, 1, 47.

[3] A. M. Yousefi, H. Atsbha, Polym. Eng. Sci. 2009, 49, 229.

[4] S. Tanoue, T. Kajiwara, Y. Iemoto, K. Funatsu, Polym. Eng. Sci. 1998, 38, 409.

[5] S. Tanoue, Y. Iemoto, Polym. Eng. Sci. 1999, 39, 2172.

[6] S. I. Tanifuji, S. Kikuchi, J. I. Takimoto, K. Koyama, Polym. Eng. Sci. 2000, 40, 1878.

[7] J. F. Agassant, P. Avenas, J. P. Sergent, P. J. Carreau, Polymer Processing: Principles and Modeling, Hanser Publishers, Munich/ Vienna/New York 1991.

[8] K. Atalik, R. K. Keunings, J. Non-Newton Fluid Mech. 2002, 102, 299.

[9] B. Meulenbroek, C. Storm, A. N. Morozov, W. van Saarloos, J Non-Newton Fluid Mech. 2004, 116, 235.

[10] A. M. Grillet, A. C. B. Bogaerds, G. W. M. Peters, F. T. P. Baaijens, J. RheoL 2002, 46, 651.

[11] G. B. Andeen, J. Appl. Mech. 1970, 37, 1170.

[12] R. Fattal, R. Kupferman, J. Non-Newton Fluid Mech. 2004, 123, 281.

[13] R. Guenette, M. Fortin, J. Non-Newton Fluid Mech. 1995, 60, 27.

[14] E. O. A. Carew, P. Townsend, M. F. Webster, J. Non-Newton Fluid Mech. 1993, 50, 253.

[15] A. Di Mascio, R. Broglia, R. Muscari, Comput. Fluids 2007, 36, 868.

[16] M. J. Crochet, A. R. Davies, K. Walter, Numerical Simulation of Non-Newtonian Flow, Elsevier, Amsterdam 1984.

[17] R. Keunings, Fundamentals of Computer Modeling for Polymer Processing, Carl Hanser Verlag, Munich 1989, p. 402.

[18] R. I. Tanner, N. Phan-Thien, J. Non-Newtonian Fluid Mech. 1985, 18, 143.

[19] M. A. Hulsen, R. Fattal, R. Kupferman, J. Non-Newton Fluid Mech. 2005, 127, 27.

[20] S. Osher, R. P. Fedkiw, J. Comput. Phys. 2001, 169, 463.

[21] D. Enright, R. Fedkiw, J. Ferziger, I. Mitchell, J. Comput. Phys. 2002, 183, 83.

[22] M. Sussman, P. Smereka, S. Osher, J. Comput. Phys. 1994, 114, 146.

[23] B. Yang, J. Ouyang, Q. Li, Z. Zhao, C. Liu, J. Non-Newtonian FluidMech. 2010, 165, 1275.

[24] L. Ville, L. Silva, T. Coupez, Int. J. Numer. Meth. Fluids 2011, 66, 324.

[25] A. Bonito, J. L. Guermond, S. Lee, Int. J. Numer. Methods Fluids 2016, 80, 53.

[26] A. Bonito, M. Picasso, M. Laso, J. Comput. Phys. 2006, 215, 691.

[27] K. K. Kabanemi, J. F. H6tu, Rheol. Acta 2009, 48, 801.

[28] T. J. R. Hughes, L. P. Franca, G. M. Hublert, Comput. Meth. Appl. Mech. Eng. 1989, 73, 173.

[29] M. Behr, T. E. Tezduyar, Comput. Meth. Appl. Mech. Eng. 1994, 112, 3.

[30] O. Pironneau, Finite Element Methods for Fluids, Wiley, Chichester 1989.

[31] F. Yurun, M. J. Crochet, J. Non-Newton Fluid Mech 1995, 57, 283.

[32] G. Karypis, V. Kumar, Proc. 8th SIAM Conference on Parallel Processing for Scientific Computing, Minneapolis, Minnesota, USA 1997.

[33] S. P. Timoshenko, J. M. Gere, Theory of Elastic Stability, McGraw Hill, New York, NY 1988.

[34] R. I. Tanner, J. Polym. Sci 1970, 8, 2067.

[35] M. J. Crochet, R. Keunings, J. Non-Newton Fluid Mech. 1982, 10, 339.

[36] N. Merah, F. Saghir, Z. Khan, A. Bazoune, Plast. Rubber Compos. 2006, 35, 226.

Caption: FIGURE 1 Predicted sheet wrinkling during free extrusion from a rectangular die with a large width-to-thickness ratio of the die exit cross section [chi] = 125

Caption: FIGURE 2 Sketch of the problem: fluid emerging from a rectangular die and fluid domain [[OMEGA].sub.L] in gray

Caption: FIGURE 3 A sketch of the model structure for energetic wrinkling analysis, A, flat shape of an extruded sheet, B, deformed shape after wrinkling

Caption: FIGURE 4 Predicted parison shape for a Newtonian fluid, A, global view of the extruded sheet, B, thickness swell in a vertical cross section, C, close up view of thickness swell, D, width swell in a horizontal cross section, E, close up view of the width swell

Caption: FIGURE 5 Predicted parison shape for a Newtonian fluid, A, without sag, B, with sag effect

Caption: FIGURE 6 Experimental data of steady shear viscosity [eta] of HDPE Lupolen 4261A at 200[degrees]C and its fitting curve for the PTT model, also the first normal stress difference Ni/2

Caption: FIGURE 7 Predicted parison shape as a function of extrusion time at Weissenberg number Wi = 40 and for a die exit cross-sectional aspect ratio [chi] = 12.5, using a viscoelastic PTT model: front view of the extruded sheet with sag effect

Caption: FIGURE 8 Predicted parison shape at Weissenberg number Wi = 40 and for a die exit cross-sectional aspect ratio [chi] = 12.5, using a viscoelastic PTT model, A, perspective view of the extruded sheet with sag effect, B, cross-sectional shape of the extruded sheet at a distance of 60 mm from the die exit showing thickness and width swelling

Caption: FIGURE 9 Variation in the first normal stress difference [N.sub.1] = [[tau].sub.zz] - [[tau].sub.yy] along the central line, in the flow direction, A, Wi = 40 and die exit cross-sectional aspect ratio [chi] = 12.5, B, Wi = 27.8 and die exit cross-sectional aspect ratio [chi] = 62.5, C, Wi = 29 and die exit crosssectional aspect ratio [chi] = 62.5, D, Wi = 40 and die exit cross-sectional aspect ratio [chi] = 62.5

Caption: FIGURE 10 Predicted parison shape at time = 160 seconds, Weissenberg number Wi = 27.8 and for a die exit cross-sectional aspect ratio [chi] = 62.5, using a viscoelastic PTT model, A, extrudate shape, B, cross-sectional shape of the extruded sheet at a distance of 60 mm from the die exit showing thickness and width swelling

Caption: FIGURE 11 Predicted parison shape at time = 160 seconds, Weissenberg number Wi = 29 and for a die exit cross-sectional aspect ratio [chi] =-62.5, using a viscoelastic PTT model, A extrudate shape, B, cross-sectional shape of the extruded sheet at a distance of 60 mm from the die exit showing sinusoidal wrinkle deformation of small amplitude, C, close up view of the central region highlighting sinusoidal wrinkle deformation

Caption: FIGURE 12 Predicted parison shape at time = 160 seconds, Weissenberg number Wi = 40 and for a die exit cross-sectional aspect ratio [chi] = 62.5, using a viscoelastic PTT model, A, extrudate shape, B, cross-sectional shape of the extruded sheet at a distance of 60 mm from the die exit showing sinusoidal wrinkle deformation of large amplitude

Caption: FIGURE 13 Predicted parison shape as a function of extrusion time at Weissenberg number Wi = 40 and for a die exit cross-sectional aspect ratio [chi] = 62.5, using a viscoelastic PTT model: front view of the extruded sheet showing sinusoidal wrinkle deformation

Caption: FIGURE 14 Predicted parison shape as a function of extrusion time at Weissenberg number Wi = 40 and for a die exit cross-sectional aspect ratio [chi] = 62.5, using a viscoelastic PTT model: perspective view of the extruded sheet showing sinusoidal wrinkle deformation

Caption: FIGURE 15 Predicted parison shape at two different Weissenberg numbers, for a die exit cross-sectional aspect ratio [chi] = 62.5, using a viscoelastic PTT model, A, extrudate shape at Weissenberg number Wi = 40, B, cross-sectional shape of the extruded sheet at a distance of 60 mm from the die exit showing sine-like wrinkles of large amplitude at Weissenberg number Wi = 40, C, extrudate shape at Weissenberg number Wi = 52, D, cross-sectional shape of the extruded sheet at a distance of 60 mm from the die exit showing cosine-like wrinkles of large amplitude at Weissenberg number Wi = 52

TABLE 1 Values of Phan-Thien-Tanner parameter model for HDPE Lupolen 4261A at 200[degrees]C PTT model parameters Values [epsilon] 0.015 Relaxation time [lambda] at 200[degrees] C (s) 17.3 Zero shear viscosity [eta] at 200[degrees] C (Pa.s) 1.29*[10.sup.5]