Achraf Kallel

Elie Hachem

Francesca Rapetti

Yves Demay

Jean-François F Agassant

Stability analysis of a polymer film casting problem

Keywords: cast-film, tracking strategy, capturing strategy, draw resonance, domain decomposition technique

teaching and research institutions in France or abroad, or from public or private research centers.

INTRODUCTION

The cast-film process is commonly used to produce polymer films or coatings, especially for packaging applications. It consists in extruding a molten polymer in a flat die (gap 2h 0 ∼ 1 mm, extrusion velocity u 0 ∼ 1cm/s; width W ∼ 1m) and then stretching the polymer sheet by a chill roll whose linear velocity can reach u f = 1m/s, to produce a polymer film of thickness 2h f (of the order of 10 µm). The stretching distance L between the die and the roll is generally of the order of several cm (see Figure 1) [START_REF] Agassant | Mise en forme des polymères[END_REF]. When the hot polymer film (temperature around 200 • C) is in contact with the cold roll (temperature around 20 • C), it is cooled down rapidly due to its very thin thickness in such a way that stretching is really imposed on the distance L.

The key parameters of this cast film process are the Draw ratio Dr = u f /u 0 and the shape factor A = L/2h 0 . Due to the high viscosity of molten polymers (η = 10 3 P a.s) the final gap of the die 2h 0 needs to be sufficiently large to limit the pressure required to deliver the appropriate extrusion Figure 1. Typical cast film process and leading parameters visualization velocity u 0 . It is thus necessary to impose a high Draw ratio to produce the target film thickness 2h f = 10µm and this results in several processing problems.

First, the film width will shrink during stretching which limits the film thickness reduction (h f /h 0 < 1/Dr). In addition, this "neck-in" phenomenon will result in an inhomogeneous final film thickness even if the initial sheet thickness at die exit is homogeneous. The resulting localized overthicknesses on each side of the film (commonly called "dog-bone" phenomenon) make it necessary to trim the lateral borders of the film and this will again reduce the width of the produced film. Both phenomena are amplified when the Draw ratio or/and the stretching distance increase. One can try to reduce the stretching distance but this induces complex technological problems and generally the shape factor is higher than 10 (L > 1cm).

Moreover, periodic instabilities on both the film width and the film thickness are observed as soon as the Draw ratio is higher than a critical value (Dr > Dr c ) (see Figure 2). Figure 2. Typical width and thickness periodic instabilities in cast film processing [START_REF] Barq | Instability phenomena in film casting process[END_REF] These defects represent obvious limitations for the productivity of the process and a lot of scientific and technological effort has been devoted to limit the "necking" and "dog-bone" phenomena and to postpone the development of periodic instabilities. Increasing the drawing distance L and applying an appropriate cooling air jet along this drawing distance allow to increase the critical Draw ratio Dr c but at the same time both "necking" and "dog-bone" phenomena are enhanced. One can try to reduce these defects by reducing the stretching distance and it was speculated that the critical Draw ratio would remain constant or even decrease (which means that periodic instabilities would appear at a lower stretching velocity) but surprisingly, the opposite was experimentally observed. This paper aims at understanding these observations.

STATE OF THE ART

Due to the industrial importance of the cast film process, a large range of models exists (see for example [START_REF] Demay | An overview of molten polymer drawing instabilities[END_REF]) and the goal here is not to give an exhaustive description for each of them but to underline the main common features connected with this study. The different models deal principally with one or several of the following difficulties:

• behaviour of the molten polymer (Newtonian or viscoelastic),

• description of "necking" and "dog-bone" phenomena,

• determination of the critical Draw ratio,

• description of the transition from a shear flow to a stretching flow.

The first model is a simple but interesting 1D membrane model [START_REF] Pearson | Mechanical principles of polymer melt processing[END_REF]. It assumes a constant width and a Newtonian behaviour. If variables x and t denote respectively stretching direction and time, the unknown functions h(t, x) and u(t, x) denote respectively the half thickness of the web and the mean value of the velocity in the xdirection over the thickness. The steady model has an analytical solution which obviously does not account for "necking" and "dog-bone" phenomena. Stability of the steady flow is studied using linear stability analysis. It consists in calculating the dominant eigenvalue (two complex conjugates eigenvalues) of the stationary solution. This Newtonian isothermal constant width situation leads to the well known critical Draw ratio Dr c = 20.2 below which all eigenvalues have a negative real part which means that any flow perturbation will decrease and finally vanish as a function of time [START_REF] Fisher | Finite amplitude stability and draw resonance in isothermal melt spinning[END_REF]. Let us mention that for Dr > Dr c , a Hopf bifurcation occurs: the steady flow becomes unstable when a time-periodic solution appears. This elegant and rapid stability analysis has been then applied to more realistic and complex situations (viscoelastic constitutive equations, non isothermal flow) in the case of fiber spinning which is basically governed by the same system of equations (see [START_REF] Demay | An overview of molten polymer drawing instabilities[END_REF] for a detailed bibliography). Of particular interest, this approach points out that increasing the drawing distance and introducing film cooling will stabilize the process (leading to a higher critical Draw ratio) and this is consistent with experiments. The drawback of this method is that computing the time dependent solution is tedious and possible only in a small neighborhood of the critical Draw ratio. Let us point out that, as the transition from a shear flow in the flat die to a stretching flow between the die and the chill roll is neglected, the critical Draw ratio is overestimated: this is easily understood having in mind that swelling (in the y direction) occuring at the die exit is neglected.

The following are membrane models using mean values of velocity components u(t, x, z) and v(t, x, z) in the x and z directions respectively (stretching and shrinking directions) as well as thickness e(t, x, z) ([6]; [START_REF] Silagy | Stationary and stability analysis of the film casting process[END_REF]; [START_REF] Debbaut | Viscoelastic effects in film casting[END_REF]; [START_REF] Kim | Transient solutions of the dynamics of film casting process using a 2d viscoelastic model[END_REF]; [START_REF] Shin | Transient and steady-state solutions of 2d viscoelastic nonisothermal simulation model of film casting process via finite element method[END_REF]). The important contribution of these 2D membrane models is a description of "necking" and "dog-bone" phenomena. They were first developed based on a Newtonian behaviour but they were then extended to account for viscoelastic fluids. In this case, the stability of the steady flow is determined using transient simulations and a convenient time-space numerical scheme. Growth or decay of a small perturbation of the steady flow can be observed. It was found that the critical Draw ratio depends on the initial film width 2W 0 and the stretching distance L (see Figure 1).

A difficult point is to take into account the transition from a shear flow in the flat die to a stretching flow between die and chill roll (Figure 3). This can be done for large initial width by assuming a constant width of the web and considering a 2D Stokes free boundary problem in the (x, y) plane of Figure 1 ( [START_REF] Souli | Finite-element study of the draw resonance instability[END_REF]; [START_REF] Fortin | Numerical simulation of coextrusion and film casting[END_REF]). For instance, since the thickness of the polymer film and the stretching distance are negligeable comparing to the width, it is possible to assume that the flow is invariant according to z (the position in the transverse direction). Generalization to a viscoelastic fluid remains challenging. In this paper, two complementary methods for the 2D problem in the (x, y) plane will be proposed and compared, first for a long stretching distance (A ≫ 1) to check the validity of the results with the existing membrane models, and then for a short distance (A ∼ 1) which corresponds to a new development.

• A front-tracking method is first developed. A stationary solution is computed in the polymer flow domain (Figure 3). The interface with air is determined to satisfy the kinematic interface steady equation. The linear stability method is then applied to this interface location in order to predict onset of Draw resonance. • For the direct numerical simulation, both polymer flow and the surrounding air are considered and the time dependent interface, between polymer and air, is captured by solving a suitable transport equation for a level-set function.

The front-tracking method coupled with a linear stability analysis is able to determine the onset of Draw resonance at a low computational cost for various processing conditions. Scilab programming language is used and it results in a computation code that runs in few minutes on a laptop computer.

FRONT-TRACKING METHOD

The two dimensional stretching flow of a (highly viscous) Newtonian fluid is considered. As this flow presents an unknown interface with air, the classical and well posed Stokes problem is coupled with a kinetic interface equation allowing to determine the air-polymer interface at each time step. This interface becomes unstable at high stretching rates. Let us point out that as the Reynolds number is small (Re ≪ 1), this instability is not inertial. The development of instability at high values of the Draw ratio is due to the nonlinear terms of the kinematic interface equation. The finite element computation of the steady flow is first described. The linear stability is studied in section 4 by introducing the time dependency in the kinetic interface equation. The comparison of linear stability results to a time dependent simulation using a capturing method will be detailed in section 5.

Governing equations for the steady flow

The considered fields are velocity U , rate of strain tensor ǫ, stress tensor σ and pressure p:

U = u v ; ǫ = 1 2 ∇U + ∇U t ; σ = σ xx σ xy σ xy σ yy
As the objective is to describe the transition from a purely shearing flow in the flat die to a purely stretching flow (Figure 3), the considered geometry includes the final part of the flat die and the stretching area till take up. The change of boundary conditions at die exit induces a geometric singularity which is responsible for computational difficulties. Let us note 2h 0 the gap of the flat die. Extrusion occurs at x = 0. The flow in the die is considered on a length h 0 (from x = -h 0 to x = 0) and the contact of the fluid with the wall is sticky. After extrusion the liquid is stretched by imposing at x = L a take-up velocity u L largely greater than the mean extrusion velocity u 0 . The flow domain Ω and its boundaries are presented on Figure 3. Let us precise the boundary conditions:

• The velocity components vanish on the upper wall Γ 1 of the flat die.

• A plane Poiseuille flow is imposed on Γ 2 .

• On the symmetry axis Γ 3 , the normal component of the velocity vanishes (v = 0) as well as the shear component of the stress tensor (σ xy = 0). • On Γ 4 (take up), the x-component of the velocity is imposed (u = u L ) and the shear components of the stress tensor vanishes (σ xy = 0). • The interface with air H is described as the graph y = h(x). On this part, the interface balance equation writes:

σ - → n = 0 .
The stress tensor σ for a Newtonian behaviour reads (η p is the polymer viscosity) :

σ = 2η p ǫ -p Id (1)
Due to the large value of η p , gravity and inertia are neglected, and the balance equation becomes :

∇ • σ = 0 (2)
Finally, as the fluid is incompressible, mass conservation holds:

∇ • U = T r ǫ = 0 (3)
Equations ( 1)-( 3) form the classical Stokes problem which is well posed with the previously described boundary conditions on ∂Ω for a given air-polymer interface (i.e for a given function y = h(x)). A finite element strategy is used to compute velocity components u, v and pressure p (see paragraph 3.2). Let us point out that as inertia is neglected, the considered equations are unchanged if the time dependency is introduced. At this step, a supplementary condition is needed to determine the steady interface H. As the non-miscibility condition writes for a steady flow U • - → n = 0 [START_REF] Souli | Finite-element study of the draw resonance instability[END_REF], we have:

u h ∂h ∂x = v h with u h (x) = u(x, h(x)) and v h (x) = v(x, h(x)).
The non-miscibility condition leads to a highly non linear equation for h, namely :

F (h) = -u h ∂h ∂x + v h = 0 (4)
An iterative process is used to solve this equation. First, the interface position is initialized by the analytical solution of the 1D membrane model h(x) = h 0 Dr -x L [START_REF] Yeow | On the stability of extending films: a model for the film casting process[END_REF]. Then, several fixedpoint iterations allow to make the interface position closer to the solution and thus ensure an easy convergence. Finally, the steady-state interface is obtained by successive Newton-Raphson's iterations.

The finite elements strategy

For a given air-polymer interface, the Stokes problem is classically written as a mixed problem:

Find U ∈ V and p ∈ Q such that:

a(U, V ) + b(V, p) = 0 ∀ V ∈ V (5a) b(U, q) = 0 ∀ q ∈ Q (5b) with a(U, V ) = 2 Ω η p ǫ (U ) : ǫ (V ) and b(U, p) = - Ω p div (U )
The spaces V h and Q h of test functions are convenient subspaces of Sobolev spaces V ⊂ H 1 (Ω)

2
and Q ⊂ L 2 (Ω) [START_REF] Ern | Éléments finis: théorie, applications, mise en oeuvre[END_REF]. The finite element approximation of problem ( 5) leads to a linear system for nodal values of U and p of the following form:

A B ⊤ B 0 U P = F U 0 (6)
In this matrix equation, U and P are the unknown vectors of nodal values for the velocity and pressure whereas the right hand side part F U of the equation contains Dirichlet conditions for velocity. Due to Korn inequality, the bilinear form a(., .) is coercive and hence matrix A is positive definite. However, the linear system ( 6) is not invertible for all choices of the finite element pairs for interpolating velocity and pressure. The finite element pair P + 1 -P 1 (Mini-element, [START_REF] Arnold | A stable finite element for the Stokes equations[END_REF])which satisfies the Babushka-Brezzi-Ladysenskaia (BBL) condition is used here for velocity components and pressure. According to this choice, every triangular element K is divided in three sub-triangles having as a common vertex the element's center of gravity and the discrete velocity space is enriched by a bubble function defined as the test function associated to the center of gravity of the element. This bubble function is continuous, linear on each sub-triangle and vanishes on the boundary of K. The supplementary velocity component degrees of freedom introduced by the bubble function are removed thanks to the so-called static condensation strategy which consists in writing the discrete equations of (5a) with the bubble as test function on each velocity component. Assembling on the three sub-triangles allows to compute nodal values of velocity at the center of each element as a function of nodal values of velocity and pressure at the vertices of the element. It thus leads to the following linear system :

A B ⊤ B C U P = F U F p (7)
In the following sections, A, B, C, F U and F p will be simply denoted by A, B, C, F U and F p .

Flow domain, meshing and domain decomposition technique

The flow domain between the extrusion die and the take-up is restricted by the position of the interface H which is defined by a piecewise linear interpolation of function h :

h i = h(x i ) with x 0 = 0 and x N = L
where h 0 is half of the flat die thickness and the set

h 1 • • • h N is unknown.
Since it is required to determine the onset of Draw resonance, for various processing conditions, relatively quickly (few minutes) and at a low computational cost, a direct and simple strategy is considered. For example, structured meshes are used to easily control node numbers at the air-polymer interface. Moreover, a reference domain

Ω ref = [-h 0 , L] × [0, h 0 ]
is defined and a meshing of domain Ω is then obtained by applying the following transformation, between die exit and take-up, to a meshing of Ω ref :

(x, y) -→ x, y h(x) h 0 for x ≥ 0
In addition, mesh refinement at die exit is required to deal with the pressure singularity resulting from the change of boundary conditions at die exit (sticky contact at the wall Γ 1 in the die and vanishing normal stress on H).As the effect of this geometrical singularity is well known for the so called stick-slip flow (in this case velocity and stress tensor behave respectivily as r 1 2 Φ(θ) and r -1 2 Ψ(θ)) and in order to avoid the use of an adaptative method, the mesh of the reference domain Ω ref is radially refined in the neighborhood of extrusion die exit (x = 0 and y = h 0 ). Figure [START_REF] Pearson | Mechanical principles of polymer melt processing[END_REF] shows the mesh refinement of the resulting flow domain. The suitability of this mesh refinement was checked since it allowed to capture the 19% die swelling when no draw ratio is imposed [START_REF] Mitsoulis | Simulation of extrudate swell from long slit and capillary dies[END_REF]. Obviously, in the flow region located immediately after die exit (about half a thickness after die exit), shear components which were dominant inside the die, vanish progressively whereas elongation components become largely dominant along the stretching path. In the remaining part of the flow, the velocity component v and its gradient become negligeable with respect to the velocity component u and its gradient respectively. Therefore, in this part of the domain, it is required to refine the mesh only in the x direction. Thus a domain decomposition approach, allowing to adapt the meshes to the features of the solution in the flow domain, is adopted [START_REF] Cockburn | Discontinuous galerkin methods for convection dominated problems[END_REF]. Indeed, it allows to obtain a convenient solution at low computational cost. Classical approaches when dealing with nonconforming finite approximations are the mortar element method ( [START_REF] Bernardi | A new nonconforming approach to domain decomposition: the mortar element method[END_REF]; [START_REF] Lions | The mixed mortar finite element method for the incompressible stokes problem: Convergence analysis[END_REF]) and the discontinuous Galerkin method. See also [START_REF] Canuto | Spectral methods : Evolution to complex geometries and applications to fluid dynamics[END_REF] and the references therein for more details on discretization methods for fluid problems and their theoretical analysis in complex geometries. The scope here is mostly confined to a very simple coupling approach which neither introduces additional difficulties (such as complicated function spaces or new variational terms depending on parameters to set) nor alters the significant features of the solution (such as stability, flow rate distribution) even if the coupling is not optimal with respect to the consistency error. The rule of mesh-conformity (i.e., two elements of the mesh sharing either nothing, or a vertex or a whole edge) constrains mesh refinement to propagate from the die exit singularity up to the end of the domain.

To by-pass this difficulty, two meshes that do not match on the connection line γ are used as shown in Figure 4. The coupling technique used to link components of U and P associated to nodes on γ for the two meshes is described in the following. The flow domain Ω is decomposed into two non-overlapping sub-domains

Ω k , k = 1, 2: Ω = Ω 1 ∪ Ω 2 , Ω 1 ∩ Ω 2 = ∅, Ω 1 ∩ Ω 2 = γ, (8)
such that Ω 1 contains the die exit singularity and Ω 2 the final part of the polymer film. The connection line γ is vertical and located after die exit, for example at x γ = h 0 where neither the pressure nor the velocity present too strong variations in the y direction. Then two different families of triangulations T 1 and T 2 are used on Ω 1 and Ω 2 respectively. Let us denote by U k and P k the vectors of nodal values for U and P on mesh T k , k = 1, 2. After the assembling on each local mesh, one obtains matrices A k , B k and C k . Let us point out that, since the nodes of T 1 and T 2 do not in general coincide on γ, the continuity conditions

U 1 (x) = U 2 (x) ; p 1 (x) = p 2 (x) ; for x ∈ γ (9)
have to be restored at the discrete level. This is done by interpolation. For the final linear algebraic system, unknown values on T 1 and T 2 are split. We set U = [U 1 , U 2 ] ⊤ and P = [P 1 , P 2 ] ⊤ with

U 1 = [U 1 I , U 1 γ ] ⊤ , U 2 = [U 2 γ , U 2 I ] ⊤ , P 1 = [P 1 I , P 1 γ ] ⊤ , P 2 = [P 2 γ , P 2 I ] ⊤ ,
where U k I (resp. P k I ) denotes the vector of unknown values for velocity (resp. pressure) associated to nodes x j ∈ T k with x j / ∈ γ and U k γ (resp. P k γ ) the one of unknowns for velocity (resp. pressure) associated to nodes x j ∈ T k ∩ γ. With the chosen ordering, we have :

A 1 = A 1 II A 1 Iγ A 1 γI A 1 γγ
, and

A 2 = A 2 γγ A 2 γI A 2 Iγ A 2

II

(similarly for B k and C k ). The right-hand side vectors F k U and F k p have the same block structure as the one of U k and P k . For the discrete continuity condition on γ, it is decided to consider Ω 1 as the (master) subdomain that passes the information to the (slave) subdomain Ω 2 at γ. This choice is motivated by the fact that T 1 is finer than T 2 on γ in the y-direction (see figure 4). Let us first describe this master-slave strategy for the pressure (a scalar field). The continuity condition ( 9) is imposed by interpolation on γ and it induces that the coefficients p 2 i , i = 1, m 2 γ , associated to nodes (x γ , y 2 i ) ∈ T 2 ∩ γ verify :

p 2 i = j=1,m 1 γ p 1 j ϕ 1 j (x γ , y 2 i ), i = 1, m 2 γ .
where ϕ 1 j is the piecewise linear Finite Element function associated to node (x γ , y 1 j ) on mesh T 1 . Here m 1 γ and m 2 γ are respectively the number of nodes on γ for T 1 and T 2 . In matrix form, we thus obtain

P 2 γ = QP 1 γ , with Q ij = ϕ 1 j (x γ , y 2 i ), 1 ≤ i ≤ m 2 γ , 1 ≤ j ≤ m 1 γ . ( 10
)
This can be written as follows:

          P 1 I P 1 γ P 2 γ P 2 I           = Q       P 1 I P 1 γ P 2 I       with Q =    Id 0 0 0 Id 0 0 Q 0 0 0 Id   
If the block matrix C has the structure

C =   C 1 0 0 C 2   =     C 1 II C 1 Iγ 0 0 C 1 γI C 1 γγ 0 0 0 0 C 2 γγ C 2 γI 0 0 C 2 Iγ C 2 II     (11)
according to the set of unknown [P 1 I , P 1 γ , P 2 γ , P 2 I ] ⊤ , the matrix C r for the reduced set of unknowns [P 1 I , P 1 γ , P 2 I ] ⊤ is then given by:

C r = Q⊤ C Q (12)
Let us note that if nodes of T 1 and T 2 coincide on γ, equation ( 12) achieves the assembly of triangles of T 1 and T 2 having nodes on γ.

The same strategy is applied for nodal values of the velocity components using matrices Q 2 and Q2 . They are built as follows:

U 2 γ = Q 2 U 1 γ and hence           U 1 I U 1 γ U 2 γ U 2 I           = Q2       U 1 I U 1 γ U 2 I       with Q2 =    Id 0 0 0 Id 0 0 Q 2 0 0 0 Id   
Matrices à and B are defined similarly to equation ( 11) by:

à =   A 1 0 0 A 2   and B =   B 1 0 0 B 2  
Then, the reduced matrices A r and B r and the reduced vectors F r,U and F r,p are defined by:

A r = QT 2 Ã Q2 , B r = Q⊤ B Q2 , F U,r = Q⊤ 2 F U F p,r = Q⊤ F p
With this notation, the final reduced system for the solution of the Stokes problem writes:

A r B ⊤ r B r C r U r P r = F U,r F p,r (13)
where

U r = [U 1 I , U 1 γ , U 2 I ] ⊤ , P r = [P 1 I , P 1 γ , P 2 I ] ⊤
The domain decomposition technique secures a good compromise between reducing the number of elements in the final film thickness and preserving a structured mesh. Thus, it helps optimizing the mesh by refining it locally only in the flow regions or in the directions where refinement is needed. Moreover, this technique improves the condition number of the global matrix compared to the classical mono-domain technique.

Numerical computation of the interface steady-state solution

The interface-tracking is used in the computation of flow problems with unknown interface. It requires meshes that "track" the interface and are updated at each iteration [START_REF] Tezduyar | Interface-tracking and interface-capturing techniques for finite element computation of moving boundaries and interfaces[END_REF]. To this end, the interface is represented by several connected marker points h 1 • • • h N across the stretching direction. At steady-state, these marker points must satisfy the kinematic interface equation (4). For a given set of parameters h 1 • • • h N , a flow domain is defined and meshed using the domain decomposition technique previously mentioned. Stokes problem is then solved in the flow domain using the mixed formulation and velocity and pressure fields are computed according to the boundary conditions previously mentioned. This allows to compute the discrete function F defined by:

F (h 1 , • • • , h N ) =    f 1 (h 1 , • • • , h N ) . . . f N (h 1 , • • • , h N )    with f i (h 1 , • • • , h N ) = -u h (x i , h i ) ∂h ∂x (x i ) + v h (x i , h i )
The set of N unknowns representing the steady interface are computed through solving the following set of N equations :

F (h 1 , • • • , h N ) = 0 (14)
This is achieved using an iterative process which is initialized according to the analytical 1D model solution mentioned previously. To improve this initialization, several iterations allow to get closer to the shape of the steady interface. These iterations consist in computing h by solving the following 1-D equation:

- d dx dh dx = - d dx u x, h(x) v with h(0) = h 0 and h(L) = h 0 Dr
Finally, equation ( 14) is solved numerically by successive Newton-Raphson iterations. At the beginning of each iteration, the interface is updated according to the computed velocity field and the gradient of the function F is numerically computed using the approximation (for a small enough value of δh):

∂f i ∂h j (h 1 , • • • , h N ) = 1 δh f i (h 1 , • • • , h j + δh, • • • , h N ) -f i (h 1 , • • • , h j , • • • , h N )
This strategy is relatively precise but, unlike the capturing strategy, it is unable to describe transient evolution of the interface.

In the following the steady interface is denoted by y = h(x).

LINEAR STABILITY ANALYSIS OF THE INTERFACE STEADY-STATE SOLUTION

Linearized equations

Once the interface steady-state solution is computed, the linear stability method is used to investigate its stability which allows determining the onset of Draw resonance instability. When considering an unsteady flow, the kinematic interface equation ( 4) becomes:

∂h ∂t = F (h) (15)
Since the only equation including time-derivatives is the dynamic interface equation [START_REF] Arnold | A stable finite element for the Stokes equations[END_REF], stability of the steady state h is studied by introducing a small perturbation as follows :

h(t, x) = h(x) + e Λt h(x) ( 16
)
where Λ is an eigenvalue and h(x) an eigenfunction. Inserting this expression into the dynamic interface equation and neglecting second order terms, gives the following linearized equation:

∇ h F (h) h = Λ h (17)
Thus, the stability problem is reduced to a classical eigenvalue problem as described in section 2. The gradient of function F is computed numerically using a similar strategy to the one considered in the case of Newton-Raphson iterations and the N eigenvalues are then computed using a classical Scilab routine. The steady-state interface is stable only if the real part of the leading eigenvalues (two complex conjugate eigenvalues) is negative. As a test, numerical results were compared to those of the classical membrane model by imposing a slipping condition on Γ 1 and considering a large stretching distance and the classical critical Draw ratio Dr c = 20.2 was obtained.

Case study :

L 2h 0 = 2 and Dr = 18

In this section, some results are shown for a specific case ( L 2h 0 = 2 and Dr = 18 ). Figure 5 shows a scaled image for the steady-state interface. Figures 6 and7 show the velocity iso-values in the stretching direction and the pressure distribution. These figures prove that the membrane assumption becomes valid at the final part of the film. This finding justifies the usefulness of block-structured meshes. Besides, the die exit singularity is revealed in Figure 7.

Linear stability method has been used to investigate the stability of the computed steady-state interface. The obtained leading eigenvalues are Λ = -1.89 ± i 11.1 and the steady-state interface is thus stable in this case.

DIRECT NUMERICAL SIMULATION

Extended domain and interface capturing strategy

The polymer flow domain Ω is extended to the air flow domain in the neighborhood of the airpolymer interface. The whole domain Ω e is represented on Figure 8. The subdomain Ω e \Ω is occupied by a fictitious Newtonian fluid of very low viscosity η a . The viscosity η a can be chosen as the air viscosity but the key point is that η a /η p ≪ 1 and the balance equation ( 2) ensures that the normal stress vanishes on the air-polymer interface H. Boundary conditions for the Stokes problem are extended to the new domain as follows (Figure 8) :

• On Γ 4 , the same boundary condition is extended to the air phase : the x-component of the velocity is imposed (u = u L ) • On Γ 5 , a free-surface condition is considered which means that normal stress vanishes:

σ -→ n = 0

• On Γ 6 , a non slip condition is considered meaning that velocity components vanish (u = 0 and v = 0) The Stokes problem (1)-( 3) is now solved on Ω e and combined with a convection equation for a level-set function (see [START_REF] Hachem | Stabilised Finite Element for high Reynolds number, LES and free surface flow problems[END_REF]). If φ(x, t) and E denote the distance function to H and a truncation thickness respectively, the level-set function α is defined as follows:

α =                    2E π for φ > E 2E π sin π 2E φ for |φ| < E - 2E π for φ < -E (18)
Function α satisfies (see [START_REF] Sussman | A level set method for computing solutions to incompressible two phase ow[END_REF]) the convected level-set equation :

     ∂α ∂t + u.∇α + λ s(α) ∇α -1 - π 2E α 2 = 0 α(t = 0, x) = α 0 (x) (19)
where the sign function s(α) is defined as follows:

s(α) =          1 if α > 0 0 if α = 0 -1 if α < 0 (20)
In this equation λ is a coupling constant depending on time and spatial discretisations (typically λ ≃ h K /∆t where h K is the mesh size and ∆t the time step). Following [START_REF] Ville | Convected level set method for the numerical simulation of fluid buckling[END_REF] and setting:

a = s(α) ∇α ∇α and g(α) = 1 - π 2E α
a rearranged form of ( 19) leads to the following simple convection equation:

     ∂α ∂t + (u + λa) • ∇α = λs(α)g(α) α(t = 0, x) = α 0 (x) (21)
Compared to [START_REF] Osher | Level set methods: An overview and some recent results[END_REF] and [START_REF] Sethian | Level Set Methods and Fast Marching Methods : Evoloving Interfaces in Computational Geometry[END_REF], this method is specific since we intend to restrict all the levelsets to the interface, thus avoiding different reinitialization of the distance function all over the domain and reducing computational effort. Moreover, this particular filtering of the levelset function allows to embed the reinitialization steps automatically in the transport equation. At this point the Stokes problem (1)-( 3) is coupled with equation ( 21) through the expression of viscosity. This is classically done using a linear interpolation for the dynamic viscosity as follows:

η = H(α)η p + (1 -H(α))η a ( 22
)
where H is a smoothed Heaviside function given by:

H(α) =          1 if α > ε 1 2 1 + α ε + 1 π sin πα ε if |α| ≤ ε 0 if α < -ε (23)
Here ε is a small parameter (the interface thickness) chosen such that, if h I is the mesh size in the normal direction to the interface, ε ∼ h I . The use of linear interpolation for the dynamic viscosity is advised for an interface orthogonal to the flow. However this can lead to inaccurate results for an interface parallel to the flow. In order to handle the abrupt changes at the interface, the following harmonic mean formulation has been used:

1 η = H(α) η p + 1 -H(α) η a (24)

Numerical strategy

The Stokes problem (1)-( 3) is solved using the mixed finite elements strategy (P + 1 -P 1 ) described in section (3). As equation ( 21) is an advection equation it is necessary to use upwind scheme. A finite element formulation based on the classical SUPG (Streamline upwind Petrov-Galerkin) method able to control the spurious oscillations of the advection dominated regime is used (see [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF]). The finite element formulation of equation ( 21) can be written as follows:

Find α h ∈ V h such that ∀w h ∈ W h Ω ∂α h ∂t + (u h + λa h ) • ∇α h -λs(α h )g(α h ) w h + n el e=1 Ω e τ SU P G ∂α h ∂t + (u h + λa h ) • ∇α h -λs(α h )g(α h ) (u h + λa h ) • ∇w h = 0 (25)
where V h and W h are standard test and weight finite element spaces for the P 1 interpolation. The classical Galerkin terms are represented by the first integral whereas the element-wise summation, tuned by the stablization parameter τ SU P G , represents the SUPG term needed to control the convection in the streamline direction. More detail about the use of stabilized finite element method for the convection equation and the evaluation of this parameter on anisotorpic meshes can be found in [START_REF] Coupez | Solution of high-reynolds incompressible flow with stabilized finite element and adaptive anisotropic meshing[END_REF].

Finally, the last step consists in deriving the anisotropic adapted mesh that describes very accurately the interface between the polymer and the surrounding air. Indeed, anisotropic mesh adaptation allows to capture scale heterogeneities in particular when simulating multiphase flow problems : the discontinuities and high gradients of the solution are captured with good accuracy for a reasonable number of elements. This challenging construction of the anisotropic mesh adaptation can be conveyed to that of constructing an appropriate mesh tensor by means of a discretization error analysis. All details of the algorithm can be found in [START_REF] Coupez | Solution of high-reynolds incompressible flow with stabilized finite element and adaptive anisotropic meshing[END_REF].

Therefore, we start by applying an a posteriori error estimator based on the length distribution tensor approach and the associated edge based error estimator to ensure an accurate capture of the discontinuities and an oscillation free solution. The key idea of the approach is to produce extremely stretched and highly directional elements under the constraint of a fixed number of nodes. All boundaries and evolving inner interfaces are then captured automatically in particular for high viscosity ratio. Figure 9 shows the modified level-set function and the distribution of the viscosity along the polymer-air interface. Snapshots of the adapted mesh at different time steps are depicted in Figure 10. As expected, the algorithm progressively detects and refines the mesh at the moving polymer-air interface leading to a well respected shape in terms of curvature, angles, etc. The zoom at the take-up position illustrates the sharp capture of the interface and the right orientation and deformation of the mesh elements (longest edges parallel to the boundary). This yields a great reduction of the number of triangles and consequently a reduction in the computational cost. Moreover, the anisotropic mesh adaptation takes into account only the levelest function since the background mesh remains isotropic. Remeshing is carried out every 10 time steps and the solution is linearly interpolated on the new mesh. It is worth mentioning that the large width of the interface is related to the interface thickness and also the use of a size gradation to ensure a smooth transition, which is needed in the case of monolithic formulations.

COMPARISON OF THE TWO NUMERICAL METHODS

The two proposed numerical methods are compared in a specific case : A = L 2h 0 = 5 and Dr = 10.

These strategies are compared on the basis of the obtained steady-state interfaces and also with respect to stability. Figure 11-a compares the shape of the interface at the steady-state : both methods are able to converge to the same solution. This is also confirmed by comparing the velocity profiles along the x-axis (Figure 11-b). The interface stability is investigated using the linear stability method and direct simulation.

In the case of tracking strategy, linear stability analysis results in the following leading eigenvalues Λ = -1.15 ± i 10.39 which means that the process is stable.

In the case of capturing strategy, the transient stability analysis is performed by introducing a small perturbation to the steady-state solution. The transient response due to this perturbation is monitored by the direct simulation, taking into account both domains, the polymer fluid and the surrounding air. The evolution of the final film thickness h(L) is depicted in Figure 12. Convergence to a steady value of the final thickness is observed within a precision of 1% and this shows that the numerical method is reasonably conservative. The real part of the leading eigenvalues is estimated from the exponential envelope of the transient response while its imaginary part is estimated from the time period of oscillations. Leading eigenvalues Λ = -1.18 ± i 9.88 are obtained for the same conditions which is very close to the leading eigenvalues computed with the tracking strategy.

This finding shows that both strategies are able to give the same stability results using two different stability analysis techniques.

NUMERICAL EXPERIMENTS

Since the tracking strategy is faster and easier to implement, it is used in the following to investigate the process stability.

Steady-state interface evolution with the shape factor

A = L 2h 0
Figure 13 shows the shape of the steady state interface obtained for a Draw ratio of 18 and two different shape factors. In the case of A = 15, a 4% die swelling is observed. When a planar Poiseuille flow is considered and no Draw ratio is imposed, a Newtonian die swell of about 19% is numerically established [START_REF] Mitsoulis | Simulation of extrudate swell from long slit and capillary dies[END_REF]. However, by increasing the Draw ratio, die swell decreases. In the case of A = 5 which means a stretching distance three times shorter, no die swelling is observed. Therefore, die swelling is much influenced by the stretching distance and it vanishes for low values of A.

Interface stability as a function of the shape factor

A = L 2h 0
A Draw ratio of 18 is considered and the interface stability is investigated for three different shape factors. Table I shows the dominant eigenvalues Λ obtained in each situation. The case A = 2 is very stable because the real part of its dominant eigenvalues is negative and relatively far from zero. Besides, the case A = 5 is stable but it is just at the limit of stability. Thus, reducing the stretching distance (at a fixed Draw ratio) improves significantly the stability of the process. This result confirms experimental studies (performed at Bostik company) showing that when the stretching distance is of the same order of magnitude as the die gap, the process becomes more stable at very high Draw ratios.

In the case A = 15 the real part of the dominant eigenvalue is close to zero. This means that the critical Draw ratio is close to 18. Since its shape factor is relatively high, this result can be compared to the critical Draw ratio Dr c = 20.2 given by a constant width membrane model. The difference between the two solutions is not surprising since a different set of boundary conditions is used in each model. In the membrane model, a plug flow is considered at die exit which does not account for the flow inside the die. Nevertheless, the two critical Draw ratios have the same order of magnitude. Table I. The dominant eigenvalues for different shape factors

Interface stability as a function of the Draw ratio

Three different Draw ratios are considered and interface stability is investigated as a function of the shape factor A = L 2h 0 . Figure 14 shows the evolution of the real part of the dominant eigenvalues Λ as a function of the shape factor A. For all Draw ratios, the real part of the dominant eigenvalues increases with the shape factor A. For instance, for a low Draw ratio Dr = 18, the real part of the dominant eigenvalues remains negative until important stretching distances : It becomes slightly positive at A = 15. However, for a higher Draw ratios (Dr = 36 and Dr = 100), the real part of the dominant eigenvalues becomes positive for a shape factor A of the order of 3. By reducing the shape factor to A = 2, the real part of the dominant eigenvalues becomes significantly negative and of the same order of magnitude for all Draw ratios which means that the process becomes very stable regardless of the Draw ratio. Interface tracking method and direct numerical simulation give similar results for the stationary interface of a Newtonian film stretching process. The tracking strategy coupled with a domain decomposition meshing technique makes it possible to apply an efficient linear stability method to determine the stability of the process. The direct numerical capturing method gives similar results which validates the tracking method. Nevertheless it necessitates important computation resources and it is thus difficult to use it for a systematic investigation of the extrusion and stretching conditions. This has been done with the tracking method to investigate cast film processing conditions. For large stretching distances, the real part of a pair of complex conjugate eigenvalues becomes positive for a critical Draw ratio whose value is very close to the critical Draw ratio of the membrane model. When reducing the stretching distance, the real part of the dominant eigenvalues decreases which corresponds to a marked stabilizing effect and this has been observed experimentally, on a pilot cast film line at Bostik company, by reducing the distance between the extrusion die and the chill roll. For further developments more realistic viscoelastic constitutive equations need to be considered.

Figure 3 .

 3 Figure 3. Geometry of the 2D flow domain (only one half of the flow domain is considered)

Figure 4 .

 4 Figure 4. Structured mesh transformation

Figure 5 . 0 = 2 Figure 6 .Figure 7 .

 50267 Figure 5. Steady-state interface ( L 2h 0 = 2 and Dr = 18 )

Figure 8 .

 8 Figure 8. Geometry of the extended 2D flow domain

Figure 9 .

 9 Figure 9. (a) The level-set function, (b) The mixed viscosity

Figure 10 .

 10 Figure 10. Anisotropic adapted meshes at different time steps

Figure 11 .

 11 Figure 11. (a) The shape of the interface, (b) The velocity profile along the x-axis

Figure 12 .

 12 Figure 12. Transient stability analysis

Figure 13 .

 13 Figure 13. The shape of the interface for different shape factors

  ± i11.94 -0.27 ± i12.03 -1.89 ± i11.1

Figure 14 .

 14 Figure 14. Interface stability for different Draw ratios and different shape factors

ACKNOWLEDGMENT

The authors would like to thank BOSTIK SA and TOTAL GROUP for technical and financial support.