Preface

The objective of the present work is the analysis of coextrusion processes by numerical simulation based on phase-field modeling of stratified confined flows. The study of such flows is motivated by the presence of complex phenomena appearing in a vast range of industrial operational coextrusion conditions due to the differences in the components properties and their viscoelastic behavior. The basic idea in coextrusion is to combine several layers of different polymers in a common die, to form a unique product with enhanced properties. However, the existence of fluid stratification in the die is responsible of a severe distortion of the interface between the fluid components, causing a loss of efficiency for the whole process. Experimental data show that, even if a stratified initial configuration is imposed at the die entry, one fluid eventually encapsulates the other in most of the flow condition analyzed. The intrinsically three-dimensional nature of this phenomenon has required the development of a three-dimensional flow solver based on the finite volume discretization of the Navier-Stokes equations for incompressible and isothermal flow, together with differential nonlinear constitutive equations (Giesekus, PTT models). The presence of two fluid phases is taken into account by a phase field model that implies the solution of an additional scalar equation to describe the evolution of the interface on a fixed Eulerian grid. This model, unlike others of the same family, has a thermodynamic derivation and can be physically interpreted. The proposed method is tested against experimental data and solutions already available in literature and a study of coextrusion in rectangular dies is performed to identify the dependence of encapsulation on the flow parameters L'ensemble des travaux présentés dans cette thèse porte sur la simulation numérique des procédés de coextrusion par un modèle d'écoulement stratifié basé sur la méthode du champ de phase. L'avantage technologique offert par la coextrusion réside dans la possibilité de combiner des matériaux ayant des propriétés physiques très spécifiques dans un produit unique. Toutefois, les différences rhéologiques entre les divers matériaux sont elles mêmes responsables d'un phénomène de distorsion de l'interface séparant deux couches adjacents. Les données expérimentales en coextrusion bicouches montrent que, en raison des différences de viscosité et d'élasticité entre le deux composants, le fluide moins visqueux encapsule l'autre et le passage d'une configuration stratifiée à une encapsulée comporte une perte de qualité du produit final. Ce phénomène, dit d'enrobage, représente donc un sujet de très grande actualité dans la recherche industrielle et la compréhension des mécanismes le générant sera utile pour l'amélioration des procédés de mise en forme des polymères. La nature intrinsèquement tridimensionnelle de l'enrobage a requis le développement d'un code pour la simulation en trois dimensions basée sur la méthode des volumes finis pour la discrétisation des équations de Navier-Stokes pour les écoulement incompressibles et isothermes couplées avec une loi constitutive différentielle non linéaire (modèles de Giesekus ou PTT). La présence de deux fluides est prise en compte par une équation scalaire supplémentaire décrivant l'évolution de l'interface sur un maillage fixe. Cette équation offre une interprétation physique précise car elle est dérivée de la thermodynamique de séparation de phase d'un fluide binaire.

Le modèle proposé est validé par confrontation avec les résultats expérimentaux et numériques disponibles dans la littérature. Une étude numérique de la coextrusion en filière rectangulaire est effectuée afin de mettre en évidence les facteurs influençant l'enrobage et la nature de son origine.

de l'enseignement supérieur (CINES) for granting the allocation of computing resources; Gilmar Mompéan and Laurent Thais from the Laboratoire de Mécanique de Lille for providing the initial version of the numerical code upon which this work is based; Bruno Vergnes from CEMEF for providing the experimental data and Alexandre Delache for the numerous and fruitful discussions on the matter of fluid mechanics.

I would like to personally thank the members of the Examining Committee for the time they spent in reviewing my work and for attending the Thesis defence. Among them a particular attention goes to my supervisors Benoît Blottière, Jacques Guillet and Emmanuel Leriche, whose precious support throughout the course of my Ph.D. has been a cause of professional but also individual growth.

My days in France would not have been easy at all if I had not had great people by my side. All of you who I have met during what I will probably remember as the three coldest years of my life, will have a special place within my heart. Nonetheless I feel like I have to mention Taline Boyajian, Riccardo Puragliesi, Roula Fares, Ausilia Campanaro and Imène Boussetouan for being there in the good days and the not-so-good ones. Thank you! You guys were great with me and I will always be grateful for meeting you! A huge "thank you" also goes to the people from my laboratory. Your support till the end meant a lot and the memories of the moments we shared together will always stay with me.

Finally, I have to say "thank you" to my family and friends in my hometown for cheering me up when I was having a hard time. The greatest love is the one that surpasses every physical distance. Aspect ratio = 3 (c). Aspect ratio = 4 (d). Mauffrey, 1998 [39] . . . The self-crimping characteristics of the produced fiber are function of the distri- bution of the two components in the molten state within the spinneret hole. It is evident that in order to control the amplitudes and the frequencies of the crimps of the solidified fiber, the interfacial configuration must be controllable as well. For this reason the textile industry had focused its attention in the experimental study of bicomponent coextrusion in circular dies since 1970.

Contents

Flow direction Inflow

Multilayer coextrusion can also be employed to produce thin flat films or tubes.

Each layer contributes a specific property like oxygen or moisture barrier, heat sealability, stiffness or flexibility, chemical inertness, electrical resistivity as well as optical properties like color or transparency. Tubes can be further processed by blow molding to produce bottles or tanks while films and plates can be manufactured into several shapes by thermoforming (cups, trays, blisters...) or serve as coating layers in overmolding processes (automotive industry).

Coextruded films are more economical than the ones produced by the conventional lamination process and the former technique also allows much thinner layers. Moreover, cast film coextrusion presents the advantage over film-blowing in that the film is drawn down and formed onto a rotating chill roll allowing a faster cooling and thus a higher production rate.

There exists two kind of coextrusion lines that are widely employed by polymer processing industry: multi-manifold dies and single-manifold dies. The firsts allow using polymers with liberally different properties since the layers come into contact with each other only after exiting the die. At the same time, when more than three layers are needed, multi-manifold dies are difficult to design and expensive to operate. The logical alternative is represented by a single-manifold die where up to several hundreds of layers can be extruded simultaneously. The drawback is that the interaction between the layers can give rise to interface irregularities entailing poor product quality. The rest of this chapter is dedicated to the review of principal aspects concerning experimental studies as well as to an overview of the main theoretical results and numerical methods that have been adopted so far for the simulation of encapsulation in coextrusion flows.

Review of the experimental studies

The experiments on coextrusion reported in literature are divided in two distinct groups based on the historical period to which they belonged. In the years between

1972-1976: The circular die

The first observation of the phenomenon of viscoelastic encapsulation ever reported in literature was presented in the work of Southern and Ballman [START_REF] Southern | Stratified bicomponent flow of polymer melts in a tube[END_REF] at the American Chemical Society Polymer Division Meeting in 1972 and published later in 1973.

In this paper three causes are identified as possible reasons for the layer rearrangement observed in coextrusion: the elasticity difference between the components; the preferential wetting of the tube wall by one component and the viscosity difference.

Other factors, reported later in other works, include geometric quantities like the flow rate ratio; the initial repartition, i.e. layer thickness distribution; the wall surface roughness and the cross-section shape.

Most of the papers published between 1972 and 1976 indicate the viscosity mismatch

on the contact interface between two different polymers as the main reason for encapsulation, even though the simultaneous development of a second order fluid theory, discussed here in section (1.2), seemed to point out that the mechanism of the interface deformation is governed by the normal stress differences which are non-zero in polymer melts. The first and second normal stress differences, hereafter denoted as N 1 and N 2 , are defined with respect to a local coordinate system formed by the velocity direction (1), the velocity gradient direction (2), and the third orthogonal direction (3). In this system:

N 1 = σ 11 -σ 22 N 2 = σ 22 -σ 33 (1.1)
Where σ is the Cauchy stress tensor. The difficulty to match the experimental observations to the existing theory was probably due to the lack of accuracy in measuring N 1 and N 2 . The source of error is due to the fact that the measure of N 2 involves the difference of two large numbers in order to obtain a small number and thus the result is extremely noisy. Many authors were not even able to confirm the exactitude of the sign of the measured N 2 (Lee et al. [START_REF] Biing | An Experimental Study of Rheological Properties of Polymer Melts in Laminar Shear Flow and of Interface Deformation and Its Mechanisms in Two-Phase Stratified Flow[END_REF]).

The device used to measure the encapsulation state in [START_REF] Southern | Stratified bicomponent flow of polymer melts in a tube[END_REF] is a a 0.05 inches diameter per 3 inches length capillary rheometer in which the apparent wall shear rate 4Q/πR 2 is varied from 7.45 to 149 reciprocal seconds by adjusting the global throughput Q.

A dye additive in polystyrene B permits the the definition of the interface shape and position, since in the examined laminar flow conditions the diffusion of the additive across the interface is negligible. The solidified polymer samples are extracted from the rheometer and the encapsulation is then measured as the fraction of the tube wall wet by the colored component. This procedure is established as a standard and used later on by the majority of the authors.

When extruded alone, the swelling of the two commercial polystyrenes A and B, measured as the ratio between the extrudate diameter and the capillary diameter, is approximately the same, with values between 1.2 and 1.5. This quantity is taken as a measure of the melt elasticity. For the shear conditions analyzed the die swell ratio is always close to unity, which is why the elasticity difference is not considered as a driving factor in encapsulation in this case. Preferential wetting, i.e. surface tension, is also ruled out because the capillary effects between the two polystyrenes White, 1974 [36] fluid possessing the strongest (negative) N 2 should display a convex interface, while the experiments show the exact opposite. Based on this observation, the authors concluded that the normal stress differences only produce lower order effects and confirmed the results of the previous papers in stating that the phenomenon is basically viscosity-driven. Nevertheless, no plausible explanation is given about the mechanism determining the interface motion.

In the experiments of Han [START_REF] Han | A study of bicomponent coextrusion of molten polymers[END_REF] the apparatus consists of a feeding block used to merge the two melt streams into a common die. This solution offers the advantage of providing a constant and controllable flow rate for each component. The interface shapes observed for the solidified and cross-sectioned samples, do not suffer from the same irregularities characterizing the samples from the capillary rheometer. As previously noted by the above mentioned authors, the lower viscosity LDPE encapsulates the more viscous PS. Interestingly enough, years later, Sizaire [START_REF] Sizaire | Numerical study of free surface newtonian and viscoelastic flow[END_REF] used the material data provided in the original paper for the shear viscosity, shown in fig.

(1.7(a)), and the first normal stress difference to fit the parameters of a single-mode Giesekus model. These are then employed to determine the second normal stress difference which, in contrast with the results of Lee and White [START_REF] Biing | An Experimental Study of Rheological Properties of Polymer Melts in Laminar Shear Flow and of Interface Deformation and Its Mechanisms in Two-Phase Stratified Flow[END_REF], is higher for the PS than for LDPE for an equal shear stress. The fact that the numerical simulations of Sizaire nicely agree with the interface profiles shown by Han denotes the An additional observation that tends to invalidate the criterion of the mismatching viscosities is that for Newtonian fluids the interface remains flat, although it is shown to posses little stability [START_REF] Dharod | Experimental and theoretical study of the extrusion of two-phase molten polymer systems[END_REF]. For instance, the stratified configuration of oil and water evolves towards the more energetically stable core annular flow (CAF) where the water migrates towards the tube walls. A similar explanation is suggested for the viscoelastic fluids for which the flow would naturally tend to a lower dissipation state, however the mechanism of the CAF formation is due to the interfacial instability of the stratified configuration while for polymer melts it seems more natural to interpret the layer rearrangement as a manifestation of the viscoelasticity of the fluids.

A dual mechanism hypothesis

In 1975 Everage [START_REF] Jr | Theory of stratified bicomponent flow of polymer melts. ii. interface motion in transient flow[END_REF] advanced the hypothesis that the interface motion in coextrusion is the result of a dual mechanism. This idea comes from the observation of the extrudate bending at the exit of the die. In fact, in the stratified flow between two parallel plates, in reason of the different viscosity, the interface shifts towards the less viscous and higher shear rate region of the die, as also reported in [START_REF] Dharod | Experimental and theoretical study of the extrusion of two-phase molten polymer systems[END_REF]. On exiting from the channel the interface returns to its initial state as the velocity changes from a parabolic to a uniform profile. This is represented schematically in fig. (1.10).

Since the less viscous component has a higher mean velocity, it also experiences a The existence of a first regime is easily explained by the difference in the viscosities of the components provoking a pressure gradient imbalance acting on the interface. fluid properties are the same on both sides of it. This deformation needs very long dies to develop, since the magnitude of secondary flows is always much smaller than that of the main flow, and it is referred to as elastic encapsulation as opposed to the encapsulation due to the viscosity jump on the interface.

Some interesting results are included in the thesis of Texeira [START_REF] Teixeira | Étude expérimentale et numérique du phénomène d'enrobage en coextrusion[END_REF] and later Mauffrey [START_REF] Mauffrey | Étude numérique et expérimentale du phénomène d'enrobage dans les écoulements de coextrusion[END_REF] who performed a series of experiments involving bicomponent coextrusion in rectangular dies of variable aspect ratios. The polymer used include several kind of polystyrenes and low density polyethylenes. As a preliminary study, the influence of temperature is investigated by coextruding the same polymer in a side-by-side rectangular die in which the two melt streams come in contact at different temperatures while the wall of the die is kept at the temperature of the "cold" layer.

The high temperature layer encapsulates the low temperature one but the interface shape is significantly different from the one observed in the bicomponent isothermal case, as shown in fig. ). In fact due to the considerable viscosity difference between the two polymer melts the less viscous PA is not able to push the PE away from the corner where there is a "dead zone" in the flow. In general, after the corner point the encapsulation degree increases less rapidly.

The initial layer repartition is also investigated. The latter should respect as much as possible the flow rate ratio between the layers, otherwise the interface is subject to a brusque deviation from its initial position due to the rapid flow adaptation. The fact that this phenomenon has only influence on the flow immediately downstream the die entrance, suggests that varying the layer repartition while keeping a constant flow rate ratio would only affect the first regime. Finally, the effect of different surface roughness does not seem to affect the encapsulation evolution since the curves tend to collapse into a single one for the four different kind of surfaces examined. This result confirms that encapsulation is unrelated to the presence of a slipping velocity at the wall.

Although it is now widely accepted that the viscoelasticity of the fluids, and in particular the presence of a second normal stress difference, is responsible for recirculating flows, this has been regarded as a different phenomenon than the one occurring in encapsulation and therefore viscous encapsulation has been distinguished from elastic encapsulation. There has not been any experimental evidence supporting the fact that these two phenomena may be two manifestations of the same mechanism, Aspect ratio = 3 (c). Aspect ratio = 4 (d). Mauffrey, 1998 [START_REF] Mauffrey | Étude numérique et expérimentale du phénomène d'enrobage dans les écoulements de coextrusion[END_REF] Figure 1.15: Coextrusion degree for the polyethylene-polyamide pairing evidencing the presence of a plateau. Mauffrey, 1998 [39] except maybe for one unpublished work, a master thesis, cited in [START_REF] Sizaire | Numerical study of free surface newtonian and viscoelastic flow[END_REF]. Nonetheless, in recent years in both theoretical and numerical studies, the "flow instability hypothesis", advanced by Everage, Han and others, has been progressively abandoned in favor of the hypothesis that the interface motion is always caused by an imbalance of the second normal stress difference acting on it. However the latter point still needs to be consolidated and validated through experimental investigation.

Theoretical Analysis 1.2.1 The minimum viscous dissipation principle

The first historical interpretation to the phenomenon of encapsulation in bicomponent stratified flows is based on theoretical approach relying on the variational principle known as the minimum viscous dissipation (entropy production) principle (PMVD) also referred to as the Helmholtz-Korteweg principle [START_REF] Sih | Multiscaling in Molecular and Continuum Mechanics: Interaction of Time and Size from Macro to Nano: Application to biology[END_REF]. The first formulation of the principle, given by Helmholtz in 1868, states that:

For a steady flow in a viscous liquid, with the speeds of flow on the boundaries of the fluid being given steady, in the limit of small speeds, the currents in the liquid so distribute themselves that the dissipation of kinetic energy by friction is minimum. This implies that in the limit of Stokes flow the energy dissipation in a specified region of flow Ω is less than that in any other possible solenoidal velocity field assuming the same fixed values on the boundaries ∂Ω.

The Stokes equations for an incompressible isothermal flow can be expressed as:

∇ • τ -∇p = 0 ∇ • V = 0 (1.2)
where ∇• denotes the divergence operator, V is the velocity and p the pressure.

When paired with compatible boundary conditions on ∂Ω and appropriate constitutive laws for τ , the system (1.2) admits the same solutions of the variational problem associated with the minimization of the functional:

I = Ω (D • τ -p∇ • V ) dv (1.3)
With D being the strain rate tensor. If the analysis is restricted to the infinite flow between two parallel plates at distance h and power law fluids, the system (1.2) can be written as:

d dy η dw dy - dp dz = 0 (1.4)
Where η = K| dw dy | n-1 . Equation (1.4) can be completed by assigning a no-slip condition at the walls (i.e. w = 0), and specifying a value for the pressure gradient or alternatively the flow rate (i.e.

h/2

-h/2 w dy). The corresponding functional to be minimized is:

I = h/2 -h/2 K dw dy (n+1) + (n + 1) dp dz w dy (1.5)
When dealing with two-phase flows, the problem can be formulated mathematically by writing the system (1.4) separately for the two regions occupied by each fluid and imposing the following kinematic conditions on the free surface:

V • n = 0 [V • t 1 ] = 0 [V • t 2 ] = 0 (1.6)
The relations (1.6), express the vanishing velocity normal to the interface and the continuity of the components along the directions t 1 and t 2 tangent to the interface.

Since the interface is a material surface, it must be in equilibrium. Hence the following dynamic conditions must be satisfied at any point on the interface:

[σ • n] = 2σ t κn [σ • t 1 ] = 0 [σ • t 2 ] = 0 (1.7)
with σ t the surface tension and κ the local interface curvature. occupying the same volume fraction in the channel and found that:

I A-B-A < I A-B (1.8)
which, accordingly with the experimental observations, indicates that the sheathcore configuration with a more viscous core is the energetically favored solution.

When considering the first term of the integral, denoted as I ′ and corresponding to the viscous dissipation, one finds that:

I ′ A-B-A > I ′ A-B (1.9)
indicating that the viscous dissipation in the sheath-core configuration is greater than the one in the side-by-side configuration. This surprising result can be ex- (1.9), the following relation can be deduced:

I ′′ A-B-A < I ′′ A-B (1.10)
which, given the fact that the pressure gradient is equal (and negative) for both configurations, leads to the conclusion that Q A-B-A > Q A-B . Consequently, for pressure driven two-phase flows the PMVD should be formulated as follows:

(PMVD1). For a fixed pressure gradient, the interface configuration observed experimentally is the one corresponding to the flow that maximizes the global flow rate.

The apparent discrepancy with the original formulation is explained by the fact that the viscous dissipations can only be compared if the velocities on the boundaries are the same which means that the flow rates for each fluid must be fixed.

If that is the case, the interface condition is no longer an arbitrary parameter. Let us assume that the flow rate of each fluid and their viscosities are fixed. For instance, let us take

Q A A-B-A = Q A A-B = Q B A-B-A = Q B A-B and η B /η A = 2.
For both configurations, there exists only one free surface position satisfying the equality of the pressure gradient on both sides of the surface. Evidently, a jump of the pressure gradient across the free surface cannot exist since it would imply a pressure jump growing infinitely in the flow direction. Therefore, the solution of the velocity profile is unique and corresponds to the one minimizing the entropy production and the pressure gradient jump across the interface. In these conditions the viscous dissipation of the sheath-core configuration is less than the one of the side-by-side configuration and consequently the pressure gradient needed to sustain the flow is also inferior (see fig. (1.17)). This consideration leads to a second possible formulation of the PMVD for two-phase flows:

(PMVD2). For fixed flow rates, the interface configuration observed experimentally is the one corresponding to the flow that minimizes the viscous dissipation or, alternatively, the pressure gradient.

Note that the formulation of the principle is based on the a-priori knowledge, from the experimental observation, of what the fluid configuration is most likely to be, although it gives no explanation on the mechanism driving the layer rearrangement.

The first regime of encapsulation

The pressure balance idea, expressed in the previous section, can be used to elucidate, at least for Newtonian fluids, the effect of protrusion of the high viscosity fluid into the low viscosity one reported in many experimental studies like [START_REF] Dharod | Experimental and theoretical study of the extrusion of two-phase molten polymer systems[END_REF]. Let us consider the side-by-side symmetric configuration A-B with a static fixed interface separating the two fluids. If the flow rate ratio and the layer repartition are both 1 : 1 but the viscosity ratio is not unitary, the fluid with the higher viscosity will also need a higher pressure gradient to flow, since for Newtonian fluid the latter is proportional to the viscosity for a given flow rate and the section area. If we now let the interface free to move, this will react to the pressure imbalance by moving towards the region initially occupied by the less viscous fluid until the pressure equilibration is restored. At the equilibrium the flow is then characterized by a less viscous but higher shear rate region and a more viscous but lower shear rate region.

The same concept is extended by Everage in [START_REF] Jr | Theory of stratified bicomponent flow of polymer melts. ii. interface motion in transient flow[END_REF] to the semi-infinite case with an The following conclusion can be drawn from the observation of fig.

(1.18):

1. The interface deviation from its initial position is an increasing function of the viscosity ratio.

2. The entry length, defined arbitrarily as the length required for the free surface displacement to be within 1% of the steady state position, is also an increasing function of the viscosity ratio. However its magnitude is always of the order of one diameter.

Even if the assumption of rectilinear flow is not valid within the entry length and the analysis is carried out for Newtonian fluids, these results basically confirm the observation of the existence of a first regime of encapsulation which depends on three parameters being:

1. The viscosity ratio 2. The flow rate ratio

The initial layer repartition

A similar procedure as the one seen above is adopted by White [START_REF] White | Theory of interface distortion in stratified two-phase flow[END_REF] to calculate the cross-section interface profiles in circular sections for Newtonian fluids with different viscosity ratios. The results confirmed the fact that, starting from an initially flat interface, the more viscous component expands in the region previously occupied by the less viscous one, showing a convex surface. However, the little curvature observed in the computed interface profiles is most likely a consequence of keeping the contact point fixed at the wall. In real situations, contact lines are able to slip in order to "follow" the dynamics of the free surface motion, therefore the computed interface profile may probably overestimate the interface curvature.

More importantly it should be noted that in these results, there is nothing indicating the movement of the less viscous fluid as to encapsulate the other. This is indicative of the fact that the second regime of encapsulation is related to the viscoelasticity of the fluids and cannot be interpreted as an effect of a viscosity difference alone.

The second regime of encapsulation

Since a generalized Newtonian fluid model cannot explain the presence of a second regime of encapsulation, a fluid model capable of reasonably describe at least the normal stress effects should be used to extend the analysis. In [START_REF] Dharod | Experimental and theoretical study of the extrusion of two-phase molten polymer systems[END_REF] a second order model is used to demonstrate analytically that single-phase rectilinear flow (i.e. no secondary flow) in a straight channel can only exists if at least one of the following conditions is met:

1. The second normal stress difference coefficient of the fluid Ψ 2 = N 2 / γ2 is zero, with γ being the value of local shear rate.

2. Ψ 2 is proportional to the viscosity of the fluid 3. The derivatives in one direction perpendicular to the main flow direction are uniformly zero

The third condition corresponds either to bidimensional or axisymmetric flow. In all other circumstances, the flow will display secondary motions due to the fluid rheology and a non radially-symmetric geometric configuration. In addition, White analyzed the two-phase symmetric stratified flow in circular sections in the hypothesis of rectilinear flow, equal viscosities but unequal Ψ 2 , and found out that:

1. For an initially flat interface the normal stress jump across the interface is related to N 2 .

2. The dynamic equilibrium condition [σ nn ] = 0 fails to be satisfied at least in one point of the interface.

From the latter result it can be observed that no flat interface can therefore exist for two fluids having different second normal stress differences. According to the sign of [σ nn ] and the fact that N 2 is negative for polymer melts, the fluid having the greater absolute N 2 should present a convex surface. The argument used is similar to the one used by Tanner to explain the convexity of the free surface of a polymer stream flowing in an open channel in [START_REF] Tanner | Some Methods for Estimating the Normal Stress Functions in Viscometric Flows[END_REF].

The same idea is exploited by Khan and Han in [START_REF] Ashfaq | On the interface deformation in the stratified two-phase flow of viscoelastic fluids[END_REF], who used a perturbation method to analyze the interface deformation for the Coleman-Noll second order fluid. Under the hypothesis that the magnitude of the transverse velocities are small enough when compared to the main velocity, the authors decoupled and solved separately the momentum equation in the flow direction first, and in the transverse directions then. Appropriate kinematic conditions are considered on the flat interface. However, for a fixed interface, a closed solution can be obtained by simply imposing the continuity of the velocity and shear stress without imposing the continuity of the normal stress. After computing the solution, if the normal stress across the interface is not balanced, the component with a greater normal stress is supposed to be pushing into the one having a lower normal stress pointing into the interface.

By expanding [σ nn ] in Taylor series around the interface, the following first order approximation for the surface displacement can be obtained.

δy = -[σ nn ] ∂ [σ nn ] ∂y (1.11)
The interface profile is calculated for different viscosity ratios and elasticity ratios and the following conclusions are drawn:

1. in absence of an elasticity difference the more viscous fluid pushes into the less viscous one;

2. in absence of a viscosity difference the more elastic fluid pushes into the less elastic one;

3. in presence of both effects the viscosity plays a stronger role than elasticity in determining the equilibrium shape of the interface between the two fluids.

Khan and Han also shown that, by eliminating the pressure from σ nn , the expression (1.11) can be further manipulated into a form evidencing the dependence of δy from the jump of N 2 across the interface. Therefore, although this analysis is based on simplifying assumptions, it can be concluded that the interface deformation observable for polymer melts during the second encapsulation regime is driven by a N 2 imbalance just as the mechanism of the first regime is driven by a pressure imbalance. Since a jump of the normal stress cannot subsist at the equilibrium configuration the interface moves as to minimize [N 2 ], that is the interface is pulled towards the fluid presenting the higher normal stress.

Numerical Simulation

The numerical simulation of two-phase flows of viscoelastic fluids presents some very specific issues that make the computation of three-dimensional coextrusion flows quite a challenging problem.

The classical approach to solve two-phase flows is to use an interface-fitted mesh.

The flow in each sub-domain is determined by the solution of the stokes equations (1.2) completed by an appropriate modeling choice for the stress tensor. The coupling across the interface boundary is assured by applying the boundary conditions expressed by (1.6) and (1.7). Since the surface position is not known a priori, an iterative procedure must be adopted. Starting by a tentative initial guess for the flow field and the interface shape, which do not satisfy the kinematic and dynamic conditions, the solution is updated by using a Newton-Raphson or a Picard iteration.

Many authors have employed the so called pathlines method who decouples the free surface calculation from the calculation of the field variables through a Picard type iteration [START_REF] Karagiannis | Three-dimensional studies on bicomponent extrusion[END_REF][START_REF] Sizaire | Numerical study of free surface newtonian and viscoelastic flow[END_REF][START_REF] Takase | Three-dimensional viscoelastic numerical analysis of the encapsulation phenomena in coextrusion[END_REF][START_REF] Sunwoo | Three-dimensional viscoelastic simulation of coextrusion process: comparison with experimental data[END_REF]. From the knowledge of a tentative velocity field, the points on the interface are updated by solving:

du dx = dv dy = dw dz (1.12)
The velocity field is then recomputed and the procedure is iterated until convergence is reached within a desired tolerance. As reported in [START_REF] Karagiannis | Three-dimensional studies on bicomponent extrusion[END_REF], for this method to converge a stabilizing underrelaxation scheme must be employed with an underrelaxation parameter ranging between 0.2 and 0.8 depending on the nonlinearity of the problem.

Other methods of updating the surface position include, the spine method used in [START_REF] Karagiannis | Three-dimensional studies on bicomponent extrusion[END_REF][START_REF] Sizaire | Numerical study of free surface newtonian and viscoelastic flow[END_REF], or the line kinematic condition found in [START_REF] Sizaire | Numerical study of free surface newtonian and viscoelastic flow[END_REF].

When a transient simulation is considered, the Arbitrary Lagrangian Eulerian (ALE)

formulation is employed to compute the solution on a moving grid and in this case the kinematic boundary condition becomes:

(V -ẋ) • n = 0 (1.13)
where ẋ is the time derivative of the interface position. The mesh velocity in the ALE formulation is not entirely arbitrary because it is associated to the remeshing procedure needed due to the interface deformation. Since this procedure must guarantee that the mesh stays as smooth as possible as the solution evolves in time, the interior mesh points must also move along with the boundary points lying on the interface. The regularity of the mesh is something that is not always possible to control since for distorted interfaces the mesh elements can result as excessively stretched as to provoke a dramatical decay in the solution accuracy. This is the case when the encapsulation approaches its final state and the two contact lines join. The computational domain undergoes a change of the topology both in space and time and if the computer code is not able to continuously adapt the mesh to this topology change, the solution is stopped and the the domain must be remeshed from scratch. In [START_REF] Sizaire | Numerical study of free surface newtonian and viscoelastic flow[END_REF], the solution is actually stopped much earlier because, even for relatively small encapsulation degrees, the elements distortion in the region close to the contact line becomes excessive, as shown in fig. (1.19(a)).

Another delicate point that requires a careful handling is the movement of the contact line. In theory, the points of contact between the interface and the wall should not move if a no-slip condition is applied. In reality however, contact lines are observed to move since the no-slip condition and the same flow equations cease to be valid in the region where two fluids meet the wall. Two solutions can be adopted in order to include the contact line motion in a simulation. The first is to introduce a certain slip amount in the contact line region. This is somewhat artificial and needs an accurate tuning of the slip parameters and the a priori knowledge of what the slip direction would be, to produce realistic results. The other solution would be to extrapolate the contact point at the wall from the position of the interface in the inner domain and using the value of the static contact angle to assign the equivalent of a Neumann boundary condition to the interface.

The slip velocity and the contact angle conditions need a precise definition of the boundary surface normal direction to be applicable. This implies that both conditions become critical when the contact line encounters a discontinuity in the boundary normal, that is when there is a corner or an edge. For instance, this may cause a code crash when the problem of coextrusion in a rectangular die is considered and the contact line approaches the corner.

Surface-capturing methods can circumvent both the problems described above. The mesh is no longer deforming with the interface but the latter is resolved on a fixed Eulerian grid, as the one of fig. (1.19(b)). Each of the two phases is described by a scalar function whose mathematical form yields different methods. For the Volume of Fluid (VOF) and Phase-field method the scalar function represents or is related to the fluid concentration and the interface is treated as a moving discontinuity. In the Level-Set (LS) method the same function is given by the signed distance from the interface which is then represented by the implicit surface c (x, y, z, t) = 0, where c is the color function. The difference between the Phase-field and VOF approaches is that while the former considers a diffuse interface the latter is based on sharp surface reconstruction (as also in the LS method) from the cell averaged values of the color function. The advantage of a diffuse interface is that it can naturally take into account for wall slipping since the same mechanism controlling the interface diffusion also governs the phase diffusion at the wall, and therefore the contact line slippage. Furthermore, any topology change can be easily handled without any restriction on the connectivity of the domain. However, the model is complicated by the necessity to solve an additional scalar equation for the phase field and the fact that the latter introduces a few other parameters whose physical significance is not always straightforward. The real drawback of the method is that the presence of a diffuse interface results inevitably in a low order accuracy in the numerical method and that the convergence of the solution can become a complicated matter because it depends both on the the grid convergence, as the mesh element size tends to zero, and on the "model convergence" as the diffuse interface tends to a sharp one, that is the interface width is reduced to zero. Recently this approach has been used successfully by Yue et al. in [START_REF] Yue | Elastic encapsulation in bicomponent stratified flow of viscoelastic fluids[END_REF] to model coextrusion flow by adaptive Finite Element Method (FEM). However the flow is simulated under the simplifying assumption of quasi2d approximation, since a three-dimensional computation would be extremely costly with the fully implicit time-stepping procedure adopted in the paper. Encapsulation is in fact a slow process and takes a very long tube length to develop which implies that, in order to observe a full encapsulation, the computational domain must be considerably stretched in one direction and the simulation must be extended to very long times.

Aim of the present work

The review of the available literature on the subject of encapsulation in coextrusion The principal aim of this thesis is to validate the latter hypothesis through a threedimensional flow analysis, and to confirm the results concerning the nature and the causes of the two different regimes. The flow of polymer melts is supposed to be isothermal for the reason exposed in [START_REF] Mauffrey | Étude numérique et expérimentale du phénomène d'enrobage dans les écoulements de coextrusion[END_REF], and gravity effects are considered negligible because of the very small density differences between the polymers melts usually employed in coextrusion. Flow field solutions shall be obtained by numerical simulation. The methodology should rely on a surface capturing method and more specifically the choice of the method is turned to the Phase-field modeling approach in reason of the successful use already reported in a previous work [START_REF] Yue | Elastic encapsulation in bicomponent stratified flow of viscoelastic fluids[END_REF], and its ease of implementation. The choice of adopting Finite Volume (FV) discretization is instead dictated by historical reasons, since the present work started on the base of the work of Mompean at al. Although most parts were coded from scratch and ported to another programming language to facilitate the integration with recently developed external libraries, the code used in this thesis should be considered as a logical extension of the method described in [START_REF] Mompean | Finite volume simulation of viscoelastic flows in general orthogonal coordinates[END_REF] to the case of two-phase flow. The two principal design constraints imposed by the case in analysis in the implementation of the code, were the necessity to be able to cope with arbitrarily high viscosity jumps in the flow field and long solution times which motivated the choice of a semi-implicit time-stepping procedure for the temporal integration of the equations as opposed to the fully explicit procedure employed in the original code.

2.1

The phase-field model

Derivation of the governing equations

Phase-field methods are based on the principal assumption that the sharp interface existing between two non-miscible fluids is replaced by a thin but numerically resolvable region where the two fluid components are considered as locally miscible.

This miscibility layer is by definition a diffuse interface and its finite extension is referred to as the interface thickness ǫ . The two components can be regarded as the two phases of a binary fluid and described by means of a scalar function φ called phase-field variable. As in level-set or volume of fluid methods, φ has the meaning of a color function and it is directly related to the phase concentration of the fluid.

Cahn and Hilliard [START_REF] Cahn | Free Energy of a Nonuniform System. I. Interfacial Free Energy[END_REF] hypothesized that the temporal evolution of the phase composition is governed by diffusive dynamics and that the interface diffusive fluxes are proportional to chemical potential gradients.

The equation

∂φ ∂t = κ∆ψ (2.1)
is referred to as the Cahn-Hilliard equation (CH). The quantity ψ denotes the chemical potential while the diffusion parameter κ is named mobility. The symbol ∆

indicates the laplacian operator.

The evolution of the binary solution is driven by the minimization of a Ginzburg-Landau free energy functional Ω f dV [START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF], with f the free energy density:

f = α|∇φ| 2 + βζ (φ) (2.2)
Hence, the chemical potential form can be derived as the rate of change of the free energy functional:

ψ = ∂ ∂φ Ω f dV = -α∆φ + βζ ′ (φ) (2.3)
where the coefficients α and β are treated as phenomenological positive constants.

The first of the two contributions in eq. (2.2) corresponds to the gradient interface energy and it governs the creation of an interface of finite thickness between the two phases, while the second corresponds to the bulk energy modeling the phaseseparation and taking into account for the non-miscibility of the fluids.

Among its many applications, the CH equation can be used to describe the spinodal decomposition between two fluid phases. As depicted in figure (2.1), the evolution of an initially homogeneous bicomponent solution is characterized by phase separation and interface creation. The interface thickness ǫ and the surface tension σ t between the phases are controllable by properly choosing the appropriate parameters in the chemical potential expression. The correlation between ǫ, σ t and the coefficients α and β will be discussed in section (2.1.2).

The simplest form of bulk free energy is given by: This choice produces two stable phases corresponding to the two minima of the function (2.4) at φ = -1 2 and φ = 1 2 and a peak energy at φ = 0, as shown in figure (2.2). The concentration of each species is easily recovered from the the phase-field function as follows:

ζ = φ + 1 2 φ - 1 2 2 (2.4)
c 2 = 1 2 + φ; c 1 = 1 -c 2 = 1 2 -φ (2.5)
We shall refer to the bulk phases corresponding to the values -1 2 and 1 2 respectively as fluid 1 and fluid 2.

The equation for the phase field evolution is obtained by injecting expression (2.3) into (2.1) and adding an advection term due to the presence of a flow field V .

∂φ ∂t + V • ∇φ = -κ∇ • ∇ α∆φ -βφ 4φ 2 -1 (2.6) Equation (2.6
) is an elliptic non-linear fourth order partial differential equation describing the advected dynamics of the phase field variable φ involving phase separa- tion, interface creation and dissolution phenomena. If the velocity field is divergencefree, as it is for incompressible flows, both the advective and diffusive term can be manipulated into a conservative form. This implies that, when paired with appropriate boundary conditions, the phase-field equation ensures the mass conservation for both fluid components. For instance, if one imposes that (V φ) • n = 0 and κ ∂ψ ∂n = 0 on the boundaries, where n is the unit normal vector, the integral of φ over the domain is constant in time.

∂ ∂t

Ω φdV = 0 (2.7)
It is very important that the mass conservation property expressed by the condition (2.7) is retained by the numerical scheme in order to correctly describe the physics of two-phase incompressible flows.

To complete the model, equations for mass, momentum and energy conservation must be specified or, if the flow is incompressible and isothermal as in our case, only equation for mass conservation and momentum balance are sufficient to close the system. These two are represented by the Navier-Stokes equations.

(i) Mass conservation:

∇ • V = 0 (2.8) (ii) Momentum balance: ρ DV Dt = ∇ • σ -℘ (φ) (2.9)
In equation (2.9) the symbol D Dt denotes the lagrangian derivative, ρ is the fluid density, σ is the Cauchy stress tensor. As usual σ can be expressed as τ -pI, where τ is the stress tensor and pI a spherical part, the pressure, required for incompressibility.

The term ℘ represents the regularized surface tension. In the sharp interface description surface tension would be equal to σ t δ (xx i ), with δ the Dirac function located at the interface, while in a phase-field model it is treated as a distributed body force term function of φ.

The form of the latter function can be derived according to the following considerations:

(i) In equation (2.6) the advection term is responsible of changing amount the free energy in the mixture by thinning or thickening the diffuse interface.

(ii) The rate of change of free energy due to the action of the velocity field must be compensated by an equal rate of change of the kinetic energy, since the global energy, sum of the two, must be conserved.

The term ℘ is thus tailored in order to produce a variation of kinetic energy equal and contrary to the variation of free energy due to advection. The resulting function is

℘ = φ∇ψ (2.10)
A detailed explanation for the derivation of expression (2.10) can be found in [START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF] for compressible flows. For incompressible flows the procedure is more or less the same but an additional gradient term must be added to the momentum equation. The latter can be easily incorporated into the pressure term which, for incompressible flows, has no thermodynamical meaning but can be interpreted mathematically as a potential term whose gradient enforces the divergence-free constraint.

To recover the "true" pressure as defined for the case of one phase flow, the surface tension forcing must be recast as the divergence of a tensor Ξ whose components are given by

Ξ i,j =              α i =j ∂φ ∂x i 2 if i = j -α ∂φ ∂x i ∂φ ∂x j otherwise (2.11)
Expression (2.10) is still preferable over (2.11) because in the tensorial representation the discretization of the surface tension term becomes non trivial.

As a consequence of introducing a phase field φ, all of the fluid properties such as the density are defined as a function of the local phase concentration:

ρ = ρ 1 (0.5 -φ) + ρ 2 (0.5 + φ) (2.12)
In the present case the density differences between the two components are neglected for the reasons already exposed in section (1.4), but if the buoyancy force is to be taken into account it can be easily included by an additional Boussinesq term in the momentum balance equation in which the density is modeled as in the formula (2.12).

The same consideration holds for the global stress tensor which is determined by the linear mixing rule

τ = τ 1 (0.5 -φ) + τ 2 (0.5 + φ) (2.13)
The whole set of equation for the phase-field model reads as follow:

∇ • V = 0 ρ DV Dt = -∇p + ∇ • τ -φ∇ψ Dφ Dt = -κ∆ [α∆φ -βφ (4φ 2 -1)] (2.14)
The form of the stress tensor in each phase is to be specified by an appropriate constitutive equation. The fact of using a mixing rule for the stress tensor, instead of the material parameters themselves, allows to take into account the possibility of using completely different constitutive equations for each phase. It is worth noticing that more complex mixing rules apart from linear are also possible, like for instance exponential or power law.

Physics of the phase-field model

The most simple equilibrium solution for the system (2.14) can be obtained analytically for the case of a planar interface in absence of flow. For this case, the phase-field equation reduces to the following one-dimensional form:

α d 2 φ dx 2 -βφ 4φ 2 -1 = 0 (2.15)
A solution of (2.15) with homogeneous asymptotic Neumann boundary condition is

given by: The surface tension can be evaluated by making use of (2.16) and integrating the continuum tension forcing term in its stress form.

φ (x) = 1 2 tanh β 2α x (2.
σ t = α +∞ -∞ dφ dx 2 dx = αβ 18 (2.18)
The latter result suggests that a proper choice for the surface tension scale unit woud be:

σ t ∝ αβ (2.19)
Assuming for instance that the interface thickness is defined as the width of the region where 90% of the phase variation takes place, and using the solution (2.16), ǫ can be estimated to be 4.164 α/β. This region also contains 98.5% of the total surface tension stress.

Instead of assigning the paramters α and β, in phase-field computations it is more convenient to be able to control the surface tension and the interface thickness, because the former is a macroscopic parameter that can be easily measured for any pair of fluids by experimental techniques, and the latter can be chosen on the base of purely numerical considerations deriving from the discretization of the CH equation.

Indeed, real interfaces always have a finite thickness but since it may be smaller than the macroscopic flow scales by several order of magnitude, it is often impractical to resolve the flow with sufficient resolution to describe both scales. Therefore, in the classical formulation of the Navier-Stokes equations, interfaces and free surfaces are conveniently threated as discontinuities.

Conversely, in the phase-field approach, it is more suitable to artificially diffuse the interface so that its width would match the scale of numerical resolution. By doing so, one introduces an error in the solution deriving from the unphysical thickness of the interface but, on the other hand, the error can be easily controlled if the right behavior is given to ǫ as the numerical grid resolution parameter h tends to zero.

This consideration yields an important criterion in the choice of the ǫ parameter.

In practice, for the interface to be correctly resolved by the numerical methods, its thickness must be greater than the mesh cell size and the following constraint, based on an empirical rule [START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF], can be adopted.

ǫ h > 7 2 (2.20)
This also implies that, as h tends to zero, ǫ must also tend to zero but with slower trend, for the diffuse interface equation to converge to the sharp interface limit, while inscreasing the interface profile reasolution ǫ/h at the same time.

The coefficient κ is also significant in phase-field simulations. Like the interface thickness ǫ, the mobility parameter can have a physical relevance in the model especially for capillary flows or any flow observed at microscopic scales since it plays an important role by governing the dynamics of interfaces coalescence. For instance, if the problem of two droplets in a shear flow is considered as in [START_REF] Yu | Coalescence of droplets in viscoelastic matrix with diffuse interface under simple shear flow[END_REF],

the effect of a high mobility is to promote the interface short range interactions between the droplets, whereas for small values their coalescence is hindered. For flows at macroscopic scales, mobility ceases to play such a fundamental role and it can be conveniently tuned to reproduce the correct interface behavior of the numerical model. Since the diffusion term of the equation is non-linear, it can act in both a diffusive and anti-diffusive way in governing the dynamic relaxation of the interface towards its equilibrium profile, so it can serve as to maintain the equilibrium profile against the thickening or thinning modes due to the advection and the smoothing action of the artificial diffusivity introduced by the numerical scheme. A high mobility can over-damp the flow while a too small value will provoke an excessively smearing of the interface profile. Good values for κ can be determined empirically by numerical experimentation and depend on the problem examined. A second and more delicate condition is to be assigned for φ. Even if no slip velocity is prescribed at the wall in the NS equation, the contact line movement can still be taken into account thanks to the diffusive term in the CH equation. In fact as the interface position is perturbed at the interior of the domain the contact line moves because of the diffusion flux tangential to the wall arising from the deviation of the interface from its equilibrium profile (see figure (2.4)). This mechanism makes the contact line follow the motion of the interface in the internal region with a velocity that depends on the mobility parameter. A high value of κ means that the contact line reacts quickly to the perturbations to which it is subject while a small value implies that the movement of the contact line is lagged compared to the interface motion. L (φ) is the wall chemical potential while the value of the constant γ depends on the value of the equilibrium contact angle. For example, for γ = 0 the contact angle will be π/2. Equation (2.23) reflects a more complex condition that can be used to take into account the non-equilibrium deviation of the contact angle from its static value and the presence of a slip velocity u s as reported in [START_REF]Contact-line dynamics of a diffuse fluid interface[END_REF].

Boundary Conditions

∂φ ∂t + u s • ∇φ = D w ∂φ ∂n -γg (φ) (2.23)
It is clear that as D w grows to infinity the static contact angle condition is recovered.

A general model for the slip velocity compatible with the phase-field model is offered by the Generalized Navier Boundary Condition (GNBC) described by Quian et al. [START_REF] Tiezheng Qian | A variational approach to moving contact line hydrodynamics[END_REF].

For a Newtonian fluid, the slip velocity in the wall-tangent direction t obeys the following law:

β s (φ) u st = -η ∂v t ∂n + ∂v n ∂t + L (φ) ∂φ ∂t (2.24)
where n denotes the wall-normal direction and β s is a slipping coefficient function of the phase composition.

The interesting feature of this boundary condition is that, in the hypothesis of Stokes is replaced by the viscoelastic wall shear stress.

Although boundary conditions (2.24) and (2.23) allow a very sophisticated modeling of the physics of contact line motion, their numerical implementation is a non-trivial matter and makes them less attractive because they further enforce the coupling between the phase, the velocity and the stress fields on the domain boundary. More-over, the coefficients appearing in their mathematical expressions are not always available nor easily measurable. If the flow is to be observed at macroscopic level, boundary condition (2.22) together with NBC or even no-slip boundary condition still offer a better option.

Constitutive equations

η 0 = lim γ→0 η ( γ) η ∞ = lim γ→∞ η ( γ) (2.25)
are referred to as zero shear viscosity and infinite shear viscosity. A general model for the stress tensor is given by

τ = 2η s D + n modes i=1 τ e i f (τ e i ) + λ i ⋄ τ e i = 2η p i D (2.26)
where η s is the so called solvent viscosity, while η p is the polymer viscosity. τ e i is the polymer stress tensor relative to the i-th relaxation time λ i and D = 1 2 ∇V + ∇V T is the strain rate tensor . The operator

⋄ () indicates the Gordon-Schowalter's (GS) convected derivative ⋄ τ e i = ∂τ e i ∂t + V • ∇τ e i -∇V T • τ e i -τ e i • ∇V + ξ (D • τ e i + τ e i • D) (2.27)
In expression (2.27) the parameter ξ is an empirical coefficient measuring the extent of non-affine motion of polymer molecules in the continuum. For ξ = 1 GS derivative reduces to the upper convected derivative while for ξ = -1 the lower convected derivative is retrieved.

The zero shear viscosity predicted by this model is given by

η 0 = η s + n modes i=1 η p i (2.28)
A natural relaxation time is also defined as

λ H = n modes i=1 η p i λ i η 0 (2.29)
This two quantities are essential in the process of non-dimensionalization of the equation and their definition is independent of the choice of the function f (τ e i ) that characterizes a particular constitutive model.

Giesekus model

Giesekus model [START_REF] Giesekus | A simple constitutive equation for polymeric fluids based on the concept of configuration-dependent tensorial mobility[END_REF] is obtained with

f (τ e i ) = I + λ i a i η p i τ e i (2.30)
and ξ = -1.

In equation (2.30), I is the unitary tensor and a is a positive mobility coefficient taking into account for the non-isotropic brownian motion and hydrodynamic drag of the polymer chains in the solvent medium. Its admissible values range between 0 and 1. The presence of a solvent contribution in the expression (2.26) must not necessarily imply the presence of a physical solution, instead it is necessary to be able to predict an infinite shear viscosity asymptote especially when a single mode approximation is used. In fact the term 2η s D can be interpreted as an additional purely Newtonian "0-th" mode.

For small amplitude oscillatory shear-flow, the non linear term of the equation is negligible and the predicted loss and storage moduli are the same as for the Upper-Convected Maxwell fluid model (UCM).

G ′ (ω) = n modes i=1 η p i ω 2 λ i 1 + (ωλ i ) 2 G ′′ (ω) = η s ω + n modes i=1 η p i ω 1 + (ωλ i ) 2 (2.31)
The material functions of shear viscosity and normal stress differences coefficients for steady shear flow are derived in closed form in the original paper by Giesekus.

η ( γ) = η s + n modes i=1 η p i (1 -g i ) 2 1 + (1 -2a i ) g i (2.32) Ψ 1 ( γ) = N 1 γ2 = n modes i=1 (2η p i /λ i ) g i (1 -a i g i ) a i (1 -g i ) (2.33)
Ψ 2 ( γ) = N 2 γ2 = - n modes i=1 (η p i /λ i ) g i ( γλ i ) 2 (2.34)
In expressions (2.32), (2.33) and (2.34) the coefficient g i is defined as

g i = 1 -f i 1 + (1 -2α i ) f i (2.35)
and

f 2 i = 1 + 16a i (1 -α i ) (λ i γ) 2 -1 8a i (1 -α i ) (λ i γ) 2 (2.36)
Giesekus model predicts a non zero second normal stress difference. For low shear rates the ratio between the two normal stress differences is proportional to the mobility parameter a.

lim γ→0 N 2 N 1 = - a 2 (2.37)
When a tends to zero, N 2 also tends to zero while N 1 remains finite. Giesekus model is thus reduced to Oldroyd-b model. Furthermore, if η p is also zero, the model

Phan-Thien Tanner model

When

f (τ e i ) = exp λ i e i η p i tr (τ e i ) Exponential Form f (τ e i ) = I + λ i e i η p i tr (τ e i ) Linear Form (2.38)
The Phan-Thien Tanner (PTT) model is obtained. This model is based on the assumption that the hypothesis of affinity between the macroscopic motion of the fluid and that of the polymer chains is not always valid. For this reason the GS mixed derivative is employed without a non zero ξ parameter to simulate the slipping effect of the polymers strands with respect to the deformation of the macroscopic medium.

A further assumption takes into account the rate of creation and destruction of junctions between the polymer chains. As a consequence the polymer relaxation time and viscosity are functions of the trace of stress tensor, implying that where the flow is characterized by a high stress level, the relaxation time and viscosity decrease so that the molten polymer behaves like a Newtonian fluid. One advantage of this model is that both the elongational and non-affinity properties of the material can be controlled separately by the non-linear parameter e and ξ. On the other hand the presence of a non zero ξ introduces unphysical oscillation at the start-up of a shear flow that may affect the numerical stability of the time-stepping methods used for temporal integration of the constitutive equation. In order to identify the non-dimensional parameters governing the physics of viscoelastic two-phase flows it is necessary to proceed to the non dimensionalization of the equations. This is done by defining the characteristic units of the system considered. In this analysis all the fluid properties are scaled by the properties of the more viscous component that we indicate as the fluid 1.

Nondimensionalization of the equations

The characteristic velocity and length depend on the problem examined. For instance, in a coextrusion problem the characteristic velocity V * corresponds to the mean axial velocity while the characteristic length L * can be taken as the diameter of the section. For the rest of this section all characteristic units are indicated by the symbol () * .

The characteristic time is taken as the convective time:

t * = L * V *
In the non dimensionalization of momentum equation, the pressure is scaled by the characteristic viscous stress because for very slow flows, as in coextrusion, the diffusion term is entirely balanced by the pressure gradient. This particular choice is called Stokes scaling in opposition to Reynolds scaling for which the pressure is scaled by the dynamic pressure ρ * V * 2

p * = τ * = µ 0 1 V * L *
Characteristic surface tension and interface width are naturally defined according to the considerations developed in section (2.1.2)

σ * = αβ ǫ * = α β
By scaling each variable by the respective characteristic unit the following dimensionless groups can be made:

Reynolds number Re = ρ * V * L * µ 0 1 Capillary number Ca = V * η 0 1 σ * Deborah number De = λ H 1 V * L * Cahn number Cn = ǫ L * Péclet number P e = ǫL * V * κσ *
The Reynolds number gives a measure of the relative importance of the inertial effects with respect to the momentum diffusivity and is usually very low for flows of high viscosity fluids like polymer melts. The Capillary number measures the relative importance of the momentum diffusion with respect to the surface tension effects.

The definition given here is slightly different from the classical one since it is derived from the phase-field model. Ca can be related to the definition based on the real surface tension using the result of equation (2.18) to obtain Ca = 3 √ 2Ca.

The group λ H 1 V * L * , represents the Deborah number expressing the degree to which the fluid elasticity manifests in response to a transient deformation. For low values of De the material response is that of a viscous fluid while for high values of De the material is observed to behave like an elastic solid. When V * L * does not match the observation characteristic time of the phenomenon examined but can be interpreted as an average shear rate, the same dimensionless group is representative of the Weissemberg number W i, indicating the degree to which the nonlinear behavior of the material is exhibited. When the applied deformation is very slow, that is W i ≪ 1, the material response is approximately linear, while for the opposite situation of W i ≫ 1 the response is highly nonlinear.

The introduction of a Cahn number is made necessary by the existence of a diffuse interface and it is a measure of the spatial scale at which the phase miscibility effects take place. Real binary polymer blends have interfaces width between 2 nm and 50 nm depending on the level of compatibility [START_REF] Stamm | Interfaces Between Incompatible Polymers[END_REF][START_REF] Schubert | Influence of chain length on the interface width of an incompatible polymer blend[END_REF]. With characteristic section diameters of the order of 1 cm, the natural Cahn number is around 10 -7 . For phasefield simulations the Cahn number is chosen in such a way as to produce a interface profile that is resolvable on the current grid resolution which means that the interface width must respect the condition (2.20). Following the analysis of Jacqmin [START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF], the error in the resolution of the hyperbolic tangent profile is o (h/ǫ) n where n is the order of the scheme adopted for the numerical solution, while the error due to the introduction of a diffuse interface is o (ǫ). Equating the two errors results in a scaling relation between the Cahn number and the mesh resolution parameter hs.

Cn ∝ h L * n n + 1 (2.39)
Finally the Péclet number expresses the relative importance of the phase advection effects with respect to the phase diffusion effects. The choice of this parameter is dictated by the fact that the interface must be able to maintain its hyperbolic tangent profile against the distortion exerted by flow. This consideration implies that diffusion must balance the advection effects at this scale or equivalently that interface Péclet number must equal one.

ǫ 2 V * κσ = P eCn = 1 (2.40)
From the relation (2.40), it follows that P e ∝ 1/Cn.

Considering that the flow has two components, a viscosity ratio and an elasticity ratio can be defined.

R µ = µ 0 2 µ 0 1 R λ = λ H 2 λ H 1
Since the adopted constitutive model introduces a solvent contribution to the viscoelastic stress, a solvent fraction parameter can be defined for each component as follows:

̺ 1,2 = µ s 1,2 µ 0 1,2
The whole set of non dimensional equations reads then as follows:

∇ • V = 0 Re DV Dt = -∇p + ∇ • τ + ∇ • (2̺D) - 1 CaCn φ∇ψ Dφ Dt = - 1 P e ∆ [Cn 2 ∆φ -φ (4φ 2 -1)] τ = n modes i=1 τ e i f (τ e i ) + De i ⋄ τ e i = 2µ p i D (2.41)
In the non-dimensional form of the momentum balance equation, the stress tensor is split into its Newtonian and viscoelastic part and the parameter ̺ represents the solvent viscosity fraction defined as follows:

̺ = ̺ 1 (0.5 -φ) + ̺ 2 R µ (0.5 + φ) (2.42)
Similarly, the quantity µ p i refers to the non dimensional polymer viscosity:

µ p i = η p 1 i (0.5 -φ) + η p 2 i (0.5 + φ) /η 0 1 (2.43)
The constitutive law expressed in system of equations ( 2

De i =            De λ i λ H fluid 1 De λ i R λ λ H fluid 2 (2.44)
In addition, geometric ratios can be defined on the base of the problem analyzed.

For example in coextrusion simulations, flow rate and layer thickness ratios will also be introduced.

Résumé du chapitre 2

Le chapitre 2 présente l'ensemble des équations modélisant les écoulements biphasiques

Numerical Method

The following chapter provides the details about the numerical discretization of system of equations (2.41) in space and time. The design criteria of the code for the numerical simulation of viscoelastic two-phase flows are dictated by the presence of a large amount of unknowns. A rough estimation is given by

N unknowns = N cells × (5 + 2 × 6 × N modes ) (3.1)
where 5 are the unknowns of velocity pressure and phase and the factor 2 and 6 take into account the presence of two phases and the symmetry of the stress tensor.

For instance, if 4 modes are adopted for each component (which is common practice in viscoelastic calculations) the number of unknowns is equal to the number of cells multiplied by a factor 53.

A further critical factor is offered by the fact that the differential operators appearing in momentum balance and phase field equations are either non linear or their coefficients are time or space dependent. As a consequence, their discrete counterparts cannot be computed once and stored in memory in the setup phase of the calculation but they need to be recomputed at each time step. Matrix-free operators involve a minimum amount of update operations at each iteration but they also need an ad-hoc implementation of the matrix-vector product and matrix inversion routines. The implementation of the numerical method adopted for the solution of while all the preconditioning, function evaluation and matrix-vector multiplication routines are user-defined.

Spatial discretization

Equations (2.41) are discretized in space by standard second order finite volume approximation applied to staggered cartesian grids. With this kind of grids, staggered arrangement of the flow variables has been regarded as a better option over colocated arrangement, since it provides a better coupling between velocity and pressure.

Moreover, boundary conditions for the pressure do not need to be assigned. This particular choice is often named "Marker and Cell" (MAC) from the original scheme developed by Harlow and Welch [START_REF] Francis | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF]. The only drawback of the MAC method is in the necessity of storing four different meshes; namely, one for each velocity component, computed at the centers of their respective faces and one for the pressure which is computed at cells centers.

In the presentation of this method the subscript p is used to denote the pressure space, which is also the space used for any other scalar field, and v to indicate the velocity space. Discrete operators defining linear application mapping one space into another are denoted by a double subscript. For instance, L p,v is a linear interpolation function operating from the pressure space into the velocity space, D v,p is the discrete divergence operator applied to the velocity space and taking values in the pressure space, and so on.

Handling multiple grids is not a complicated task if the mesh has a cartesian topology but the extension to general unstructured meshes is non-trivial. The computational molecule is depicted in fig. The viscoelastic stress components are also computed at the cell centers e. A more natural choice would be to define the normal stresses at cell centers and the shear stresses at the centers of the cell edges. The space defined in such a way is denoted by the subscript s. The motivation behind this choice is that, in this way, the divergence operator D s,v , defined from the stress space to the velocity space, only involves central differences and provides an optimal coupling between the velocity and the stress fields. On the other hand in viscoelastic constitutive laws the presence of tensor products between the stress and the velocity gradient would require a great amount of interpolation operations if the stress components were computed at four different locations. Hence, the choice to store all the components at the cell centers is preferred.

As an example of application of the finite volume technique, the procedure for the discretization of the phase-field equation is illustrated. Eq. (2.6) is rewritten in a conservative form and integrated over the cell volume:

V e ∂φ ∂t dv + V e ∇ • F dv = 0 (3.2)
The term expressing the divergence of F = φV -1 P e ∇ψ is then transformed into the sum of surface integrals defined on the faces of the cell by means of the Gauss where ∆v and ∆s e±en indicate the volume of the cell and the surface of its faces.

In order to compute eq. (3.4) it is necessary to interpolate the values of φ and ∇ψ at the face-centered nodes. The numerical expression for the chemical potential at the cell centers is given by a second order finite difference approximation

ψ i,j,k = -Cn 2 [ φ i-1,j,k -2φ i,j,k + φ i+1,j,k ∆x 2 + φ i,j-1,k -2φ i,j,k + φ i,j+1,k ∆y 2 + φ i,j,k-1 -2φ i,j,k + φ i,j,k+1 ∆z 2 ] + ζ ′ (φ i,j,k ) (3.5)
The staggered representation of the chemical potential gradient at the cell faces centers is given by central differencing rule:

(∇ψ • n x ) e+e 1 = ψ i+1,j,k -ψ i,j,k ∆x (∇ψ • n y ) e+e 2 = ψ i,j+1,k -ψ i,j,k ∆y (∇ψ • n z ) e+e 3 = ψ i,j,k+1 -ψ i,j,k ∆z (3.6)
For the approximation of the advective fluxes a second order accurate total variation diminishing (TVD) scheme is employed to obtain the face values for φ:

φ f = φ c + 1 2 l (r) (φ c -φ u ) (3.7)
The quantity l (r) is a flux limiter and it is function of the parameter r = φc-φu φ d -φu , expressing the ratio of the gradients evaluated in two adjacent cells and thus representing a natural way to monitor the smoothness of the function φ. The subscripts c, u, d are defined with respect to the upwind direction as shown in fig. (3.2). A list of possible flux limiter functions is presented in the table (3.1) The flux limiter allows to switch from high order reconstruction to low order in regions where discontinuities in the solution, indicating the presence of an interface, are detected. For the scheme to be mass-conserving, the discretization of the eq.

(2.8) must be consistent with the approximation of the divergence operator adopted for the phase-field equation. A numerical approximation for the volume integral of φ can be taken as The error o (δ), depending on the accuracy of the numerical method and the tolerance criterion used to terminate the iterative solution of the mass balance and phase-field equation, can be made arbitrary small. Although controllable, this error is cumulative in time and can lead to significant loss in the mass conservation of the system for long integration times.

Surface tension

The discretization of the momentum equations is carried out in the same way and it will not be repeated. The only difference is that the three scalar equations cor-responding to the three velocity components, are integrated over cells staggered in the corresponding spatial directions. The surface tension term must be treated in a careful way in order to ensure a correct energy transfer between the momentum equation and the phase-field equation. In particular, the discrete kinetic energy variation due to the surface tension forcing must be balanced by an equal and opposite variation of discrete free energy.

The discrete kinetic energy K = 1 2 (V , V ) can be derived based on the following definition of inner product for staggered vectorial fields:

(u, v) = ∆v [ Nx i=0 Ny j=1 Nz k=1 u x i,j,k v x i,j,k + Nx i=1 Ny j=0 Nz k=1 u y i,j,k v y i,j,k + Nx i=1 Ny j=1 Nz k=0 u z i,j,k v z i,j,k ] (3.
V x i,j,k φ f x i,j,k ψ i,j,k -ψ i-1,j,k ∆x + Nx i=1 Ny j=0 Nz k=1 V y i,j,k φ f y i,j,k ψ i,j,k -ψ i,j-1,k ∆y + Nx i=1 Ny j=1 Nz k=0 V z i,j,k φ f z i,j,k ψ i,j,k -ψ i,j,k-1 ∆z ] (3.12)
In eq. (3.12) the subscripts f indicate interpolated values at the faces centers.

The discrete free energy associated to the phase field is derived from the inner product for cell-centered scalar fields. and Fortin [START_REF] Guénette | A new mixed finite element method for computing viscoelastic flows[END_REF] for the finite element methods. The contribution in the left hand side is treated implicitly to improve the stability of the method while the one on the right hand side is treated explicitly. If the temporal discretization of this two terms is consistent, their presence in the equation will not affect the global accuracy of the time stepping scheme.

Spatial coupling between velocity and stress fields can also be improved. The de- and in two dimensions. The nodes marked by the red crosses are the ones involved in the discretization of the term ∂ 2 u ∂x∂y in the u node i, j if the choice of storing the shearing components at the cells corners is adopted, while the blue ones represent the resulting stencil for the same operator in case the first derivative ∂ ∂y is computed at the cell centers (A, B, C, D) and then interpolated at the cell edges to perform the second derivative in x. It is evident in this case that the node i, j is not coupled to any of the immediate neighboring nodes in the y direction, creating an alternate lines pattern that is even more dramatic in three-dimensions. To improve the coupling, the spatial discretization is performed differently for each of the two stabilizing terms in the momentum equation in a similar way as for the Rhie-Chow [START_REF] Chow | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF] interpolation technique used to suppress the spurious pressure modes appearing in colocated FV methods. The resulting time-continuous/space-discrete Navier-Stokes equations read as follows:

D v,p V = 0 Re dV dt + N (V ) -D s,v [2 (̺ + ̺ s ) d] = -G p,v p + D s,v L p,s [τ -2̺ s L s,p d] - 1 CaCn φG p,v ψ (3.15)
The implicit term on the left hand side of eq. The spatial discretization of the constitutive law only involves the numerical approximation of the advective term appearing in the equation. Since the stress components are computed at the cell-centered nodes the same schemes adopted for the phase-field are also employed.

Temporal discretization

A second order time-stepping procedure is employed to integrate the equations in time. With reference to the general ordinary differential equation (ODE)

∂u ∂t = f (u) + g (u) (3.16)
we denote by f (u) the non-stiff part of the ODE and by g (u) the stiff one. In general the non-stiff term is also non-linear, which motivates an explicit treatment, while the stiff term must be integrated implicitly in order to avoid severe constraints on the time step size. This is not true for the phase-field equation for which the stiff term is also non linear. When applying this idea to a system of equations, the coupling terms are also treated explicitly to split the global problem in a series of small subproblems that can be solved independently.

Two numerical variable step size semi-implicit schemes are considered: a Backward Differentiation Formula (BDF):

1 ∆t n 1 + 2ω 1 + ω u n+1 -(1 + ω) u n + ω 2 1 + ω u n-1 = (1 + ω) f (u n ) -ωf u n-1 + g u n+1
and a combination of Crank-Nicolson and Adam-Bashforth (CNAB) schemes

1 ∆t n u n+1 -u n = 1 + 1 2 ω f (u n ) - 1 2 ωf u n-1 + 1 2 g u n+1 + g (u n )
Where ω = ∆t n +1 ∆t is the ratio between two successive time steps.

Both the schemes are second order accurate. CNAB scheme only requires the storage of two stages, while BDF scheme needs three. Furthermore the latter needs a special treatment for the first step since solution at time step n -1 is not available. Usually a first order Backward Euler (BE) scheme is used to start up the BDF scheme. Although less memory consuming, CNAB scheme often displays spurious oscillations.

These are produced by non-smooth initial conditions and are not damped by the numerical scheme. One possible cure would be to use BE scheme for a few steps and then switch to CNAB, but this does not always guarantee the smoothness of the solution. In the case of flows in the Stokes limit, the pressure has been observed

to oscillate between two states, differing by a constant, even when the velocity is stationary. For this reason, BDF scheme is employed for the Navier-Stokes and the constitutive equations while CNAB is used for the phase-field equation.

Treating the coupling terms in an explicit way allows to split the solution of one time step into three different stages:

(i) Update the phase field:

1 ∆t n φ n+1 - 1 2 1 P e D v,p G p,v Ψ n+1 = 1 ∆t n φ n -1 + 1 2 ω [D v,p (V φ)] n + 1 2 ω [D v,p (V φ)] n-1 + 1 2 1 P e D v,p G p,v Ψ n (3.17)
(ii) Update the stress field:

1 ∆t n 1 + 2ω 1 + ω τ n+1 e + Deτ n+1 e = (1 + ω) ∆t n τ n e - ω 2 1 + ω τ n-1 e -(1 + ω) D v,p (V τ e ) -τ e • ∇ V -(τ e • ∇ V ) T + De a 1 -β τ 2 e -2 (1 -β) d s n +ω D v,p (V τ e ) -τ e • ∇ V -(τ e • ∇ V ) T + De a 1 -β τ 2 e -2 (1 -β) d s n-1 (3.18) with ∇ V = L s,p G v,s V and d s = 1 2 ∇ V + ∇ V T .
The global stress field and solvent viscosity are found using expression (2.13) and

(2.42) with the newly computed values of φ and τ e .

(iii) Update the velocity and pressure fields:

Re ∆t n 1 + 2ω 1 + ω V n+1 + G p,v p n+1 -D s,v [2 (̺ + ̺ s ) d] n+1 = Re ∆t n (1 + ω) V n - ω 2 1 + ω V n-1 +D s,v L p,s τ e n+1 - 1 CaCn φ n+1 G p,v Ψ n+1 -(1 + ω) [ReN (V ) -D s,v (2̺ s L p,s L s,p d)] n +ω [ReN (V ) -D s,v (2̺ s L p,s L s,p d)] n-1 (3.19)
Since the advection terms are treated explicitly, the time step size is subject to a

Courant-Friedrichs-Lewy (CFL) restriction of the type max e

1 2 U e-e 1 + U e+e 1 ∆x + V e-e 2 + V e+e 2 ∆y + W e-e 3 + W e+e 3 ∆z ∆t < 1 (3.20)
In practice the time step size is found to be slightly smaller that the one imposed by the CFL condition. In fact, since the stress components evolve in time with exponential laws, the boundary and initial conditions can greatly influence the stability of the scheme. For instance an abrupt start up of the flow or discontinuities on the boundary conditions can produce overshoots in the stress magnitude that may eventually destabilize the numerical time-stepping procedure.

Iterative Solvers

Phase-Field solver

The first stage of the time-stepping procedure involves the solution of the nonlinear system of equations expressed by (3.17). For simplicity the system can be rewritten as

φ + ∆t Cn 2 P e h 4 F (φ) = rhs (3.21)
where h = min {∆x, ∆y, ∆z}. The discrete operator F (φ) is arising from the discretization of the chemical potential laplacian and is defined on the same 33-points stencil as the classical finite difference biharmonic operator. The term 1 h 2 is factorized in order to make F (φ) ≈ o(1).

The solution of (3.21) implies a Newton-Raphson iterative procedure. At each iteration the jacobian matrix must be evaluated and the resulting linearized system must be solved. This would require assembling and inverting the jacobian matrix several times per time step and can be costly in terms on computational time and memory.

In order to avoid computing the jacobian, a different approach is adopted. Instead This procedure is cycled over the unknowns until convergence is reached but instead of using it as a standalone solver, this routine is rather employed as a non linear preconditioner for the Non-linear Generalized Minimal Residual Method (NGMRES) that computes the (m + 1) th iterate by combining the previous m iterates into a minimal-residual solution after solving a small linearized optimization problem. This method is implemented in PETSc as SNESNGMRES on the base of the algorithm proposed by Osterlee and Washio [START_REF] Oosterlee | Krylov subspace acceleration of nonlinear multigrid with application to recirculating flows[END_REF] and it can be easily integrated with any problem specific preconditioning routine. In the implementation of the non linear Gauss-Seidel method proposed in this section, the unknowns are updated in a redblack pattern for better smoothing properties. It can be easily remarked that the convergence of the method depends on the value of δ = ∆tCn 2 P e h 4 since this paremeter determines the diagonal dominance of the linearized discrete diffusion operator. Indeed for δ → 0 the eigenvalues of the jacobian matrices become clustered around the unity and the NGMRES method converges in one iteration. This situation is depicted in fig. (3.4) where the convergence history of the relative residual within a tolerance of 10 -8 is shown. For values of δ less than unity one single step is enough to smooth the initial residual of a factor of 10 -14 . One would expect the convergence properties to deteriorate dramatically with increasing mesh resolution since δ scales with the fourth power of reciprocal mesh size, but since the time step is constrained to be o(h) by the CFL condition and P e and Cn must respect the conditions (2.39) and (2.40), for the phase-field to converge to the sharp interface limit, δ scales only linearly with the reciprocal mesh size. For this reason the proposed solver is found to perform very well in the range flow condition and mesh resolution considered.

Should δ be very large, the solver could be further accelerated by using non linear multigrid techniques as the Fully Approximation Storage (FAS) used by Kim et al. [START_REF] Kim | Conservative multigrid methods for cahn-hilliard fluids[END_REF].

Stokes solver

A large part of the computing time is spent in the inversion of the linear system arising from the discretization of the Navier-Stokes equations. In matrix form, and after dividing the equations by the cell volume ∆v, the system reads as follows:

   I - ∆t Re h 2 D s,v ̺G v,s ∆t Re h G p,v 1 h D v,p 0      V p   = rhs (3.24)
Since G v,s is the transposed of D s,v and G p,v is the negative transposed of D v,p , the system can be made symmetric by appropriately scaling the mass conservation discrete equation. A first manipulation of this system is made by adding the term ∆t Re h 4 G p,v ̺D v,p to the block [START_REF] Jr | Theory of stratified bicomponent flow of polymer melts. ii. interface motion in transient flow[END_REF][START_REF] Jr | Theory of stratified bicomponent flow of polymer melts. ii. interface motion in transient flow[END_REF]. This is done to cancel the mixed second derivatives in the regions of bulk flow where the viscosity is constant.

   I - ∆t Re h 2 (D s,v ̺G v,s -G p,v ̺D v,p ) ∆t Re h G p,v - ∆t Re h D v,p 0      V p   = rhs (3.25)
Since D v,p V = 0, the modified system admits the same solutions of the original one (3.24), but has better conditioning properties, and thus convergence properties, because it has a stronger diagonal dominance. System (3.25) shall be referred to as the discrete Stokes problem and it belongs to the general class of saddle-point problems.

The lack of the diagonal block (2, 2) makes it difficult to solve the system by traditional iterative methods. Splitting the Stokes problem in two subproblems for the solution of the velocity and the pressure is a common way to circumvent this difficulty. In pressure correction methods velocity and pressure are updated independently at the cost of second order error in time due to the decoupling. This error is also proportional to the reciprocal Reynolds number [START_REF] Couzy | Spectral element discretization of the unsteady navier-stokes Bibliography equations and its iterative solution on parallel computers[END_REF] and becomes significant at low Reynolds.

In general, for Stokes flows coupled methods are preferred. The majority of coupled iterative solvers rely on Krylov subspace methods with preconditioners based on the block LU factorization of the matrix in the system (3.25) and a convenient approximation of the inverse of the Schur complement (see [START_REF] Benzi | Numerical solution of saddle point problems[END_REF] for a complete survey on saddle-point problem solvers). These methods are very robust and they can be applied in a black-box fashion since they do not rely on the underlying geometry of the problem but only on the structure of the matrix. The drawback is that they require a fair amount of setup operations to build the algebraic operators used in the preconditioning routines.

In variable viscosity problems the setup time can become significantly large since all the discrete operators need to be updated at each time step. For this reason matrixfree preconditioners are adopted. Even if the problem is symmetric, the adpoted Krylov subspace methods are the Generalized Minimal RESidual (GMRES) and the Generalized Conjugate Residual (GCR) since they allow general non-symmetric preconditioning. Like for the case of the phase-field problem the routines for GMRES and GCR are based on the PETSc implementations. The preconditioner is based on geometric multigrid (GMG). The basic idea behind GMG is that since most of the iterative solvers can suppress only the high frequency components of the error, low frequencies can be effectively smoothed out by recursively interpolating the original problem on a sequence of coarser grids and applying a few iterations of a basic solver (or smoother) on each level. To achieve optimal convergence grid transfer operators and the smoother must be carefully designed. Multigrid routines are also based on PETSc implementation MG and use one V-cycle per GMRES or CGR iterate and 3 smoothing cycles per level. In the next subsections the grid transfer operators, the system scaling and the smoothers are discussed in detail.

Grid transfer operators

In this section the system (3.25) shall be rewritten as

  F G G ′ 0     V p   = rhs (3.26)
To solve (3.26) by multigrid methods, a hierarchy of grid is built by defining the number of cells N xc , N yc , N zc along the three spatial direction x,y and z on the coarsest grid, the global number of levels n l and the refinement factors r x , r y , r z . The mesh on a general level n is defined by

N xn = N xc r n x N yn = N yc r n y N zn = N zc r n z (3.27)
In this way each cell on a level n is formed by the union of r x r y r z cells on the finer level n + 1.

The action of interpolating a function from a fine level to a coarse level is called restriction and the corresponding operator is denoted by I n n+1 . The inverse operation is called prolongation and the operator mapping from a coarse level to a fine level is denoted by I n+1 n . Since the residuals r n are calculated at the same location on the MAC grid as the corresponding unknowns , there will be a total of eight different transfer operators.

Cell-centered divergence residuals d n+1 are restricted on a coarse level by taking the arithmetic average of the related values at the fine grid centers.

U n I,J,K = 1 r y r z (J+1)ry-1 j=ryJ (K+1)rz-1 k=rzK U n+1 irx,j,k (3.30)
The prolongation operator I u n+1 n is defined by nearest-neighbor interpolation along y and z, as for the cell-centered residual, and by linear interpolation along x. For the staggered node I, J, K on the coarse grid, the interpolation coefficients are calculated as c 1 = (I%r x ) /rx and c 2 = 1 -c 1 (the % symbol denoting the reminder of the integer division). The prolongated value is then given by

U n i,j,k = c 1 U n-1 i rx , j rx , k rz + c 2 U n-1 i rx +1, j rx , k rz (3.31)
The same considerations are valid for the staggered components of the residual along y and z. The global velocity transfer operators can be defined as follow:

I V n+1 n =      I u n+1 n I v n+1 n I w n+1 n      (3.32)
and

I V n n+1 =      I u n n+1 I v n n+1 I w n n+1      (3.33)
The equation for the error on a coarse grid is obtained by projecting the system (3.26) from the finer grid to the coarse and approximating the fine grid error by a prolongation of the coarse grid error.

  I V n n+1 F n+1 I V n+1 n I V n n+1 G n+1 I p n+1 n I p n n+1 G ′n+1 I V n+1 n 0     e n e n   = rhs n (3.34)
The transfer operators defined as above are perfectly consistent with the finite volume discretization of the mass balance equation, in the sense that if a velocity field is divergence-free on a coarse mesh, the same property is observed for the next finer mesh and vice versa. Furthermore the following relation holds:

G ′n I V n n+1 = I p n n+1 G ′n+1 (3.35)
The above property, together with the fact that I V = G n . This is a very important result since it means that for each level the projected divergence and gradient operators have the exact same structure as if they were discretized directly on that level. The same can not be said for the discrete diffusion operator since

I V n n+1 F n+1 I V n+1 n = F n .
This implies that the definition of F is different on each level. However, instead of projecting recursively F from the finest level to the coarsest, it is more convenient to discretize it on each level to maintain the same structure across the grid hiearchy. In order to do so, the viscosity field must also be restricted from the finest level over the whole grid sequence. The same restriction operator as the one used for the cell-centered divergence residual is used for the viscosities.

Diagonal scaling

Prior to applying the smoother at each multigrid level, the system needs to be systematically scaled in order to uniform the magnitude of all the equations in it.

Depending on the viscosity ratio of the two fluids, the diffusive term in the region occupied by one fluid can become smaller than that in the other fluid by orders of magnitude. This causes the dispersion of the eigenvalues of the Stokes discrete operator in the complex plane and consequently the loss of efficiency of the iterative solvers. To understand this, recall that according the Gershgorin theorem, the spectrum of the matrix in the linear system (3.

S p = diag{ G ′ diag S v ii -1 G ii }.
The complete scaling matrix reads as:

  S v 0 0 S p   (3.36)

Smoothers

The coupled smoothers used for the discrete Stokes problem belong to a specifically designed class for indefinite saddle problems that circumvents the problem of having zero diagonal entries corresponding to the pressure unknowns by smoothing at once a subset of unknowns corresponding to a local subdomain. This particular technique is referred to as Box-relaxation. The most natural way to form subdomains is to associate each pressure unknown to every velocity unknowns related to it. In the case of a staggered cartesian grid each cell-centered pressure unknown is coupled to six velocity unknowns corresponding to the faces of the cell considered. In this way, smoothing a whole subset of unknown is equivalent to solving a local cell based saddle-point problem. The procedure can be resumed as follows:

for i = 1 : n cells do % For each cell, define the subset of associated velocity unknowns e 1 ; e i + e 1 ; e i -e 2 ; e i + e 2 ; e i -e 3 ; e i + e 3 ] % Restrict the discrete diffusion gradient and divergence operators

s i = [e i -
F c = F (s i , s i ) G c = G (s i , e i ) G ′ c = G ′ (e i , s i ) % Compute the local residual r i = [r v (s i ) ; r p (e i )]
% Update the velocity and pressure unknowns by solving: A third candidate as coupled smoother for GMG is the Distributive Gauss Seidel method (DGS). The procedure is given by: for i = 1 : n cells do % For each cell, define the subset of associated velocity unknowns

  F c G c G ′ c 0     ∆v ∆p   =   r v (s i ) r p (e i )   (3.
s i = [e i -
F c = F (s i , s i ) G c = G (s i , e i ) G ′ c = G ′ (e i , s i ) % Restrict the discrete pressure laplacian operator G ′ G c = G ′ G s p i ,

end for

In the three analyzed cases the cells unknowns are updated in a lexicographical way.

For better smoothing properties red-black ordering can be adopted. Extension to plane or line versions for more robustness can be useful in case the unknowns are strongly coupled in some direction like, for instance, when the cells are stretched.

Parallel implementation

The PETSc Distributed Mesh object (DM) is used both to partition the threedimensional arrays holding the values for the velocity, pressure, phase and stress unknowns, and to manage the data communications between processes. The parallel solution of the Stokes problem is discussed in this section, since the latter takes the major part of the computational time spent to do one temporal iterate. The implementation of a scalable parallel multigrid solver is a very complicated matter and it is beyond the aim of this thesis. The following analysis is limited to the aspects of the parallelization encountered in our implementation of the Stokes solver without wanting to give an extensive coverage of the topic. The reader should refer to [7, ch.10] for a complete survey on more efficient parallelization techniques for multigrid.

When considering the implementation of a multigrid solver, there are many aspects that may partially limit the parallelization of the algorithm and its scalability. Domain decomposition technique is employed to divide three-dimensional arrays into smaller arrays and each is assigned to one process. The group of processes, or communicator, is given a cartesian topology and each process can be identified by three integers m, n and p.

Each subdomain needs to communicate with others through adjacent faces, edges and corners to synchronize the boundary values. These are copied in appropriate ghost cells. Since the stencil required to compute the convection term by second order TVD methods is symmetric and 5-points wide in each direction, two layers of ghost cells are required for each direction too. For a n × n × n block the number of unknowns is n 3 while the number of ghost cells is 12 × n 2 + 48 × n + 32. The schematic of a subdomain block and its ghost cells is depicted in fig. (3.6).

Since the number of operations needed for matrix multiplication, residual interpolation and smoothing is also o (n 3 ), it is clear that unless the block size is at least equal to 16, the ghost cells number is greater than the interior cells number and the communication time becomes dominant in the calculation. The situation is complicated by the fact that even if the grid is big enough to justify the parallelism, when restricted on coarser meshes for the multigrid solution the size of each subdomain becomes inevitably suboptimal and the parallelization inefficient.

Another factor limiting the parallel implementation of our multigrid solver is that the smoother acts differently when applied in parallel. In fact, in Box-relaxation the unknowns are updated in a Gauss-Seidel fashion but as the level of parallelism is increased the algorithm behaves more like the Jacobi method, since each subdomain is smoothed independently. This consideration implies that the convergence rate of the smoother decreases as the parallelism is increased.

To partially cure both problems discussed above a certain number of MPI processes can be suppressed as the grid is coarsened. In the present case, instead of suppressing the processes, the smoothing is run redundantly on subsets of processes. This is done by condensing several subdomains into one and scattering the data from one process of the newly formed subdomain to the others. In this way the domain is divided in less but bigger partitions shared by several processes. The communication time is reduced; however an additional time is spent in the scattering process and optimality can only be achieved as a trade-off between these two conflicting factors.

The domain is partitioned on the coarsest mesh and each subdomain is then refined to ensure that the same geometric partitioning is maintained across the grid hierarchy and that interpolation and restriction operations do not require communications. This choice has the drawback of producing non load-balanced blocks on the finer grids.

As an example of the solver performances, the speed-up obtained when solving the problem described in section (4.2) is estimated. The same problem is solved on three communicators consisting of 4 × 4 × 4, 6 × 6 × 6 and 8 × 8 × 8 processes and referred to as C 4 , C 6 and C 8 respectively. The speed-up factor is defined as

S i = T (G 4 ) /T (G i ),
where T (G i ) is the computational time required to solve the problem with the communicator C i .

The problem is discretized on a 256 3 grid with 5 multigrid levels. On the coarsest level the grid has 16 3 cells and the communicator size is fixed at 2 × 2 × 2, meaning that the solution is 8, 27 and 64 times redundant on C 4 , C 6 and C 8 respectively.

When the problem is solved with the communicator C 6 , the resulting coarse grid blocks can have non uniform sizes ranging from 2 3 to 3 3 because 16 is not divisible by 6. This means that on the finest grid there will be blocks of size 32 3 and 48 3 , giving a load balance ratio of 3.4 between the biggest and the smallest block. In fig.

(3.7), the speed-up for the 256 3 grid is shown in comparison with the speed-up obtained for a 384 3 grid for which the load is perfectly balanced. Scalability is better for the latter case even though in both cases the problem scales only sub-linearly for the reasons already expressed in this section.

Résumé du chapitre 3

La méthode numérique utilisée pour la solutions des équations présentées dans le A second test case is given by the stationary fully developed two-dimensional flow of an Oldroyd-B fluid between two plane plates. Oldroyd-B model can be obtained from Giesekus model, selecting a = 0. The solution for w, τ ezz and τ eyz is given by: w = 6y (y -1)

τ zy = (1 -β) γzy τ zz = 2W i γ2 zy (4.2)
where γzy = ∂w/∂y, W i = λ H γzy is the Weissemberg number and β denotes the solvent fraction as usual. The simulation are carried out on a channel geometry given by L y = 1 and L z = 16 and uniform mesh size is again taken as h = 1/2 n with n = 5, 6, 7, 8. Only half of the channel is considered for symmetry reason. For this case the phase-field is also simulated. The layer configuration is of the type A -B -A as described in section (5. shows the numerical and analytical phase-field profile at z = 16. The latter is expressed by the heaviside function H (y -h) -0.5 where h is the interface position that can be determined analytically by the procedure illustrated in section (5.3). Fig. (4.3) shows the numerical error as a function of the mesh size for the normal stress tensor profile and the phase field profile. While the scheme is confirmed to be second order accurate for the stress, the same cannot be said for φ since only only first order convergence in L 1 norm is obtained. The same result is reported by Titarev in [START_REF] Titarev | Finite-volume weno schemes for threedimensional conservation laws[END_REF] who used WENO schemes up to the eleventh order of accuracy in combination with high-order Gaussian quadrature for integration over cell faces.

It seems in fact that the presence of the discontinuity reduces the accuracy of the numerical scheme to first order no matter how high is the order of the reconstruction and integration schemes. This justifies the use of flux limited schemes that switch to first order upwind near discontinuities and also, to a certain extent, the idea of diffused interface introducing an o(ǫ) error in the model.

The rest of the chapter is organized as follows: in section (4.2) the performances of the variable viscosity Stokes solver are analyzed. The successive sections present comparisons with experimental data or existing numerical solutions for extrusion and coextrusion flows.

Analysis of the Stokes solver

In order to test the efficiency of the variable viscosity Stokes solver the SINKER benchmark test [START_REF] May | Preconditioned iterative methods for stokes flow problems arising in computational geodynamics[END_REF] is considered, consisting in a sinking high viscosity and high density cubic block in a low viscosity medium. This problem arises in mantle convection simulations where fluids are characterized by viscosity contrasts up to 10 12 .

Furuichi describes it as "one of the most difficult problems to solve by iterative techniques" [START_REF] Furuichi | Development of a stokes flow solver robust to large viscosity jumps using a schur complement approach with mixed precision arithmetic[END_REF], making it an optimal benchmark to test the robustness of a Stokes The phase field is defined by:

φ = 1 2 tanh d -0.15 √ 2Cn (4.4)
with d = max (x, y, z). In this way the steepness of the viscosity jump can be post-smoothing iterations are applied. The parameter ∆t/h 2 Re is chosen as 10 30 to simulate conditions as close as possible to the Stokes limit and the Cahn number is chosen to produce a viscosity jump over 3.5 cells since in phase-field simulation this is the minimal width necessary to correctly resolve the diffuse interface thickness.

The viscosity contrast is varied from 1 to 10 7 and results are reported in table (4 For low viscosity contrasts the three smoothers present a similar behavior. The GCR method converges even if no smoothing but only unpreconditioned GCR iterations are applied on each grid level, as a consequence of the beneficial effect of the symmetric diagonal scaling. As the viscosity ratio is increased, unpreconditioned GCR and DGS fail to converge within 50 iteration for viscosity contrasts stronger than 10 -2 and 10 -3 respectively, while FVS and DVS converge in few iterations for the whole range of ∆η examined. DVS is faster but less robust than FVS who scales better with respect to the increasing viscosity contrast. In fig. The effect of the variation of the thickness of the layer over which the viscosity jump occurs can be evinced by the results reported in table (4.2). As the layer thickness is made wider than 3.5 cells for a fixed viscosity contrast, the performances of the three smoothers are improved, while for w t /h = 2 the number of iterations required for convergence increases dramatically. This is due to the inefficiency of the grid transfer operators in interpolating a discontinuous viscosity profile over the grid hierarchy.

The discrete diffusion operator is simply re-discretized on each level on the base of the restricted viscosity field. When dealing with discontinuous coefficients a better option is represented by Galerkin-based coarse grid operators obtained from the projection of the fine grid operators. Surprisingly DVS performs better than FVS for very steep viscosity profiles.

When the Reynolds number is increased or the time step size reduced, the solver performances ameliorates considerably and even un-preconditioned GCR becomes The effect of the stretching factor is illustrated in table (4.4). For a fixed domain, consisting of a box defined by L x = L y = 1 and Lz = 10, the stretching factor is varied by increasing the number of cells along z. The considered viscosity contrast is ∆η = 10 -1 since this is a typical value for a coextrusion problem. For very stretched cells, both FVS and DVS converges very slowly. Once again FVS shows more robustness but is also more computationally expensive than DVS. Increasing the stretching factor from 1/10 to 2/5 the number of iterations needed for convergence changes from 72 and 58 for DVS and FVS respectively to 9 and 8. The improvement is so evident that the computational time also decreases despite the increasing number of unknowns. For the stretching factors 4/5 and 1 the number of iterations required is 6 for both DVS and FVS but the computational time increases because of the increased size of the problem. Note that number of iterations required for the 32 × 32 × 320 grid is the same as for the 32 × 32 × 32 suggesting a dependence on h but not on the number of unknowns.

In conclusion, the three smoothers perform similarly in the Stokes limit and for small viscosity contrasts. DVS has better global performances than the others, showing a robustness comparable to FVS but a faster execution time. For this reason it will be employed in the next sections unless stated otherwise.

Secondary flow in rectangular ducts

Secondary motions in non-circular ducts are characteristic of viscoelastic fluids exhibiting a non zero second normal stress difference. In square sections eight vortices are observed. The recirculation direction along the walls is toward the corners, while along the diagonals the flow direction is toward the section center . Increasing the width of the section with respect to height, the vortices intensity decreases and the lateral vortices becomes progressively smaller. A detailed explanation on how the phenomenology of the vortices depends on the negative N 2 can be found in [START_REF] Yue | A general criterion for viscoelastic secondary flow in pipes of noncircular cross section[END_REF],

while numerical simulations including FEM and FVM modeling are available in [START_REF] Holmes | Modeling viscoelastic secondary flows in three-dimensional noncircular ducts[END_REF][START_REF] Xue | Numerical study of secondary flows of viscoelastic fluid in straight pipes by an implicit finite volume method[END_REF][START_REF] Debbaut | On the development of secondary motions in straight channels induced by the second normal stress difference: experiments and simulations[END_REF]. In the present section simulations are carried out on a square domain corresponding to a quarter of the duct section. Symmetry is imposed on two sides while no-slip is imposed on the opposite sides. Periodicity is imposed in the third direction corresponding to the main flow direction. A schematic of the domain is (4.9(a)) and (4.9(b)) show the solution for the stress difference τ eyy -τ exx . Note that the latter does not correspond to N 2 , since normal stress differences are defined with respect to the directions 1, 2, 3 corresponding respectively to: the flow direction; the velocity gradient direction in the plane orthogonal to the flow and the third conjugate direction. At the wall, where 2 and 3 correspond to y and x (or vice versa), the second normal stress difference is equal in magnitude to τ eyy -τ exx .

Quasi-2d calculation of encapsulation

Anomalous encapsulation occurs for the symmetric three-layer configuration of the type "A-B-A". It has been observed experimentally by Khomami and Ranjbaran [START_REF] Khomami | Experimental studies of interfacial instabilities in multilayer pressure-driven flow of polymeric melts[END_REF] and reproduced qualitatively by numerical simulations in [START_REF] Yue | Elastic encapsulation in bicomponent stratified flow of viscoelastic fluids[END_REF]. The general elastically driven encapsulation criterion predicts that, in absence of a viscosity difference between the two components, the more elastic one is always encapsulated by the less elastic one. However under certain geometric conditions the inverse phenomenon may take place. If the central layer is made of the less elastic component, and its thickness is below a critical value, it will be eventually encapsulated by the more elastic one. Since elastic encapsulation is much slower than viscous encapsulation, a fully three-dimensional simulation would require a very long domain in order to be able to observe a fully developed encapsulation. For this reason, and since our fundamental aim is to test our method against an existing numerical solution for viscoelastic two-phase flow, the quasi-2d hypothesis of the original paper is retained.

A general analogy can be made between the three-dimensional spatial evolution along the main flow direction (i.e. z) and the quasi-2d temporal evolution of the interface, but the latter is not necessarily indicative of what really happens in the three-dimensional case. In two dimensions the interface deformation is due uniquely to the secondary flows in the plane of the cross-section. If the out-of-plane velocity component is sufficiently uniform across the section, the three-dimensional interface deformation spatially evolves in a similar way as the quasi-2d interface evolves in time, meaning that φ 3D (x, y, z) ≈ φ quasi2D (x, y, W t), where W is the mean streamwise velocity. Since the velocity field in the main flow direction is not uniform, reaching its maximum value at the centerline and its minimum at the wall, the action of secondary motions is greater at the wall then elsewhere. As a consequence, the real encapsulation profile will progressively deviate from the profile obtained in quasi-2d simulation as the solution advances in time, producing a completely different profiles on long times. This is the reason why quasi-2d approximation should not be used to produce quantitative results but only to understand qualitatively the driving mechanisms of the flow.

Re Ca Pe De 1,2 Cn β 1,2 a 1,2 d 2 /
The schematic for the domain is given in fig. (4.10). Again only a quarter of the section is considered for symmetry reason. The domain is discretized over 768 × 192 × 3 points, which gives a uniform mesh size slightly finer than the smallest mesh size used in the adaptive meshes of Yue. The total number of unknowns is above 1.6 × 10 6 . Solution is computed from 0 to 80 non dimensional time units with a time step size of ∆t = 0.0125. The time required is approximately 6 hours with 16 MPI processes on two processors Intel Xeon X5560, each having 4 cores and 8 threads.

The flow parameters are presented in table (4.5).

As is the previous test case, the flow is pressure driven and the flow rate is not conserved. The pressure gradient is assigned to give an initial unitary flow rate. This is done to obtain a solution that is at least qualitatively comparable to the one given by Yue even if the dynamics of pressure driven and flow rate controlled flows are essentially different. Under these conditions, the flow rate gradually adjusts to a steady state value which is approximately 5 times greater than the original value.

To interpret the results of the simulation, the characteristic numbers are rescaled by the newly computed mean streamwise velocity value and the same procedure is

Bicomponent coextrusion in rectangular dies

Numerical simulation of three-dimensional bicomponent coextrusion is performed in comparison with the experimental data obtained by Teixeira [START_REF] Teixeira | Étude expérimentale et numérique du phénomène d'enrobage en coextrusion[END_REF]. The experimental setup consists of a rectangular die with a section aspect ratio of 2 and a total length of approximately 24 diameters. The flow rates can be controlled individually for each component as well as the layers thickness at inflow section. A typical experiment is run for 30 minutes to ensure that the steady state is achieved and then inflow and outflow section of the die are blocked simultaneously by two sash blinds. Once the fluids have cooled down, the solid sample is extracted and sliced at equally distanced positions along the main axis. Since the two polymers are colored with black and white colors, the deformed interface profile can be easily read. A constant temperature of 235 ± 2 • C is maintained along the die. Even though polystyrene has temperature dependent viscosity, variations induced by temperature do not affect encapsulation enough to justify the hypothesis of non-isothermal flow.

In case the thermal properties of the components are significantly different and produce a viscosity crossover, dependence from temperature must be necessarily taken into account.

The Giesekus model parameters are reported in table (4.6) for both polymers. Knowing that the section characteristic length is 1.35 × 10 -2 m and that the flow rate is for both polymers equal to 9.28 × 10 -3 kg/s, the flow characteristic numbers are:

Re = 1.6 × 10 -5 , De 1,2 = 1.7 -0.93, Ca = 5.3 × 10 3 , R η = 0.21, β 1 = 0.

Boundary conditions

The computational domain is sketched in fig. (4.14). L x , L y , L z refer to the domain length along x, y and z, while L e is the thickness of the upper layer. The thickness of the lower layer is given by L y -L e . In our simulations L y = L x /2 = 1 and L e = 0.5.

To investigate the influence of boundary conditions at the inflow, two different cases are analyzed. In a first simulation the total length of the die L z is taken as 10 and Dirichlet boundary conditions are applied for the stress and the velocity. In particular, the velocity is given a double Poiseuille profile with a zero stress profile.

Another option is to consider the presence of an horizontal separation plate at height

Newtonian case

In a preliminary analysis, grid convergence is studied for a Newtonian equivalent system, for which the transient time and the downstream length needed for the flow to be fully developed are much less extended than for the viscoelastic case. In this way the computational effort to compute steady state solutions is considerably lowered. The test case consists of a symmetric three-layers configuration presenting a less viscous central layer. The geometry of the die is defined by L x = L y = 1 and L z = 5. Initial layer widths are 0.25, 0.5 and 0.25.

The spatial resolutions considered are 1 32 , 1 64 and 1 128 and convergence is checked through visual inspection of the interface position in the mid-plane of the die. The interface thickness is 3.5h on the coarsest mesh and scales as h 2/3 with the increasing resolution, as suggested by Jacqmin [START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF]. If the flow rate is equally partitioned among the layers and the viscosity ratio is close to unity, the resulting flow is smooth and no abrupt distortion of the interface is observed, as in fig. (4.16). Good convergence in the interface position is already achieved for h = 1 64 , despite the low order accuracy of the phase-field methods.

Conversely, if the flow is such that the interface undergoes brusque deformations, convergence is much worse. This situation occurs for instance when, in reason of a high flow rate ratio, sudden expansion of one layer into another is produced.

In fig. (4.17), the breakup of the layer in the middle is anticipated for the mesh corresponding to the lowest resolution because the interface thickness is too wide compared to the thickness of the central layer. This case shows how the convergence of the solution depends on the particular problem at hand. If the interface is too thick compared to the mesh resolution, the hyperbolic tangent profiles of two close interface may overlap and the two interface coalesce due to the short range interactions. A quantitative criterion is discussed by Beckermann in [START_REF] Sun | Sharp interface tracking using the phase-field equation[END_REF], who recommends selecting an interface width at least 4.2 times smaller than the radius of curvature of the interface.

Non-Newtonian case

In coextrusion the flow rates are balanced and the extruded materials do not present viscosity differences as strong as for the case presented in section (4.5). The adopted mesh size for all the viscoelastic simulations is 1 64 . A typical computation at this resolution, involves approximately 1.2×10 6 finite volume cells and more than 22×10 6 unknowns. The computational time to integrate the equations from t = 0 to 60 is roughly of 20 hours on the cluster SGI Altix ICE 8200 (JADE) 0 , using 64 MPI processes and with a time step size of ∆t = 0.004. The flow field computed in 0 Centre Informatique National de l'Enseignement Supérieur in Montpellier (FR) (4.20(b)), showing a good agreement with the solution that includes the presence of a separation plate. Since the spatial development of encapsulation depends on the inflow conditions, the numerical encapsulation profile at z = 5 is compared with the best matching experimental profile along the die (fourth slice in the first row), in order to prove that the interface shape is correctly predicted at least qualitatively.

For a straightforward comparison the inflow boundary condition should be carefully assigned in order to match the real inflow conditions. The aim of the present chapter is the investigation of the mechanism of encapsulation through the analysis of the flow fields obtained by fully three-dimensional simulation of coextrusion. Of the two cases considered, the two-layers coextrusion in a square section duct is examined first. The geometric symmetry of this flow allows the isolation of the effects produced by the rheological behavior of the materials from those produced by the layer distribution, since for a fixed flow rate and layer thickness the configurations A -B and B -A yield the exact same results. In addition, the square duct represents the case where the flow is most interested by secondary recirculation. The same results obtained for this case can be qualitatively extended to the case of rectangular ducts, considering that for the latter case encapsulation will be confined to the wall or, in the limit of infinite aspect ratio, it will not be present at all. For the second case, the three-layers configuration, the layer distribution plays a role as well, since, for fixed flow rates and layer thicknesses, the configurations A -B -A and B -A -B lead to different flows.

In both cases the two regimes of encapsulation are examined separately as reasonably possible. The second regime can be easily turned off if Newtonian fluids are taken in consideration, since no layer rearrangement is observed in absence of a second normal stress difference but an interface shift. Conversely, in order to observe the second regime alone, the fully developed flow corresponding to a stratified configuration is adopted as inflow boundary condition. This is equivalent to observe the flow past the region where the first regime takes place, because the pressure in the two layers is already at equilibrium at the entrance of the die. However, when using this kind of boundary condition, only the flow rates or the layer repartition can be freely chosen but not both.

It should be noted that the separation of the two regimes is not entirely possible at least when the second regime is examined and Non-Newtonian fluids are considered. This is due to the viscosity being a function of the shear rate. So when the latter is varied, that is varying Weissemberg number of the flow, the viscosity ratio assumes different values, each yielding a different flow rate ratio, even for a fixed global flow rate and initial layer repartition.

In the following sections, the flow parameters adopted are estimated on the basis of realistic values for coextrusion flows. Taking the section characteristic dimension of the order of 10 -2 m, the flow rate as 10 -3 Kg/s, the viscosity as 10 4 P a s and the surface tension as of 10 -2 N/m, yields a Reynolds number of 10 -5 and a Capillary number of 10 4 . These preliminary estimates should suffice to exclude both inertial and capillary instabilities as causes for encapsulation.

A further clarification should be made about the role of elasticity in this phenomenon. What so far has been referred to as a Deborah number, is the ratio of the fluid natural relaxation time and the flow characteristic time or the observation time of a specific phenomenon. In the process of non-dimensionalization of the equations, without any assumption on the nature of flow, this time has been taken

as the characteristic convective time, that is the ratio of the characteristic length L and the characteristic velocity W m . In pipe or channel flows the ratio W m /L is proportional to the flow average shear rate, so the dimensionless group λ H W/L is more likely descriptive of the flow's Weissemberg number, expressing the degree to which the non-linear behavior of a material is manifested in response to an imposed shear rate, rather than Deborah number. The Deborah number of encapsulation should be evaluated on the basis of the encapsulation characteristic time which, corresponds to the duration of the time-unsteady flow. This time is by far greater

Two-layers configuration

Analysis of the first regime of encapsulation

In the assumption of well balanced flow rates and layer repartition, the two-layers flow of Newtonian fluids has only one parameter to set, which is the viscosity ratio.

The latter is usually defined as being greater than one, since a high viscosity ratio naturally induces to think to a high viscosity difference. However, for the flow equations to be well scaled, it is more correct to define the viscosity ratio as being the interface is expected to become flat again. This fact is also confirmed by the observation that for two different viscosity ratios the interface shape is surprisingly similar except for the more pronounced shift induced by the smaller viscosity ratio.

The interface profiles in the symmetry plane are compared in fig. (5.5), while fig.

(5.6(a)) shows the entry length and interface position as functions of the viscosity ratio. Note that the same definition of entry length as in [START_REF] Dharod | Experimental and theoretical study of the extrusion of two-phase molten polymer systems[END_REF] is adopted. As predicted by the approximated analysis for the flow between to parallel plates, both the entry length and the interface shift are increasing functions of the reciprocal viscosity ratio. However, the interface equilibrium always occurs within a length of the order of one diameter. This result confirms that the first regime is indeed a very rapid mechanism for flows in the Stokes limit.

When considering a more general case the interface position will be the result of the combined effects of the viscosity ratio, the flow rate ratio and the layer repartition. For a specific couple of fluids any desired layer repartition can be obtained by appropriately tuning the flow rate ratio.

For each fluid the dimensionless viscous dissipation and pressure gradient are computed. The viscous dissipation is defined as

Φ = 2 Ω η (x, y, z) D : Ddv , (5.1)
where Ω = [0,

L x /2] × [0, L y ] × [0, L z ].
Both quantities are decreasing functions of the reciprocal viscosity ratio as illustrated

Simulation of interfacial instability

The fig. (5.3) also shows the evolution of the perturbation generated at the flow startup. Since the boundary conditions are applied impulsively to the initially quiescent fluids, these react with a shear wave that propagates in the flow direction.

Although the disturbance is applied in the same way for the four cases in exam, the wave evolves differently depending on the viscosity ratio. High viscosity differences seem to promote the wave deformation into a finger-like structure as the one seen in fig. (5.3(f)). This observation suggests the presence of an interfacial instability, implying the growth of small perturbations in the flow direction, giving rise to a folded surface. This kind of phenomenon is referred to as "zig-zag instability" in literature and is often encountered in coextrusion as reported, for example, by Dooley [START_REF] Dooley | Viscoelastic flow effects in multilayer polymer coextrusion[END_REF].

A complete analysis on the problem of interfacial instabilities in two-layers flows can be found in [START_REF] Yiantsios | Analysis of superposed fluids by the finite element method: Linear stability and flow development[END_REF][START_REF] Stergios | Linear stability of plane poiseuille flow of two superposed fluids[END_REF][START_REF] Yih | Instability due to viscosity stratification[END_REF]. The main results can be summarized as follows:

1. for long wavelength perturbations, reducing the thickness of the less viscous layer promotes the instability while a thin more viscous layer is favorable to the stability.

2. For short wavelength perturbations, the stability region is reduced to a tighter interval of the thickness ratio. The sole purpose of this section is to illustrate the effect producing a folded interface in order to show the difference between the latter and proper encapsulation, without investigating further on the subject of the interfacial instability.

A two-layers flow is simulated considering the following perturbation of the boundary conditions as proposed by Sizaire [START_REF] Sizaire | Numerical study of free surface newtonian and viscoelastic flow[END_REF]:

Q 1 = Q 0 + ∆Qsin (ωt) Q 2 = Q 0 -∆Qsin (ωt) (5.2)
The flow rate amplitude ∆Q = 0.35 and frequency f = 2πω = 0.25 are chosen in such a way as to produce a wave shape that resists to the damping effect introduced by the grid and the diffuse interface. For lower amplitudes or higher frequencies, the wavelength and amplitude of the resulting perturbation become of the same order of the interface width and thus the numerical error prevents the onset of the instability, because the flow cannot "sense" the perturbation unless a higher numerical resolution is used. This effect is comparable to the stabilizing effect of the surface tension. In this sense the diffused interface introduces an artificial surface tension that vanishes as Cn tends to zero and P e becomes infinite.

In order to have the same spatial and temporal resolution as in section (5.2.1), but a longer channel length (16D), only two-dimensional flow is considered. The same flow parameters are also retained. Fig. (5.7) shows the results of the simulation. When subject to an initial perturbation, the interface manifests a wavy shape. The crests of the waves are subject to a velocity gradient that tends to pull them apart in opposite directions producing pronounced folds. These become progressively thinner and more elongated until they eventually break up and form an actual intermixing layer between the two phases. In order to minimize these instabilities, Dooley suggests to

Analysis of the second regime of encapsulation

The effects of viscosity and geometry have been discussed in sections (5.2.1) and

(5. still considerably large and this, together with the unboundedness of the viscosity, makes the elongational stress large enough to become a destabilizing factor for the numerical scheme adopted for the integration of the constitutive equations. This issue could be limited to some extent by using Finitely Extensible Nonlinear Elastic (FENE) models [START_REF] Warner | Kinetic theory and rheology of dilute suspensions of infinitely extendible dumbbells[END_REF], adopting different formulations of the constitutive equations based on the matrix logarithm of the conformational tensor [START_REF] Afonso | The log-conformation tensor approach in the finite-volume method framework[END_REF][START_REF] Fattal | Constitutive laws for the matrix-logarithm of the conformation tensor[END_REF] or modeling the separation plate as having a finite thickness and a blunt edge. However, the problem is simply circumvented in this analysis by adopting a different inflow boundary condition.

For the simulations concerning the second regime, the fully developed steady flow corresponding to a fixed layer repartition is assigned as a boundary condition to the inflow section. The solution of the three-dimensional flow field is then composed of two steps:

1. The quasi2D solution for a layer repartition of 1 : 1 (meaning that the interface is kept fixed during this stage) is computed until steady state is reached.

2. This solution is imposed as inflow boundary condition as well as initial condition for the three-dimensional flow.

A schematic side view of the domain is presented in fig. (5.8(a)). The resulting flow corresponds to two streams merging with perfectly balanced flow rates. In this sense it is like observing the flow past the entry length required for the pressure equilibration. Since the entry length is always o (D) and encapsulation is usually much slower, the hypothesis of a first regime occurring without noticeable encapsulation effects is fairly reasonable. Hence, the interface is assumed as still flat at the entrance of the domain.

A further advantage of this assumption is that the viscoelastic stresses have already reached their fully developed profiles at the start of the simulation, which considerably shortens the time required to reach the steady state and also produces a faster encapsulation development as seen in section (4.5). This means that a shorter die length can be considered.

The analysis of the flow includes the following aspects:

1. Inspection of the normal stress level at the inflow section:

The normal stress difference τ xx -τ yy is considered instead of N 2 because it does not depend of the definition of a local coordinate system. As illustrated in fig. The approximation is due to the fact that v is non zero at the symmetry plane and so the flow direction is not exactly parallel to z. The normal stress difference is discontinuous at the interface since the latter is not at equilibrium and the sign of the N 2 jump across the interface is to be related to the direction of the interface displacement.

2. Analysis of the secondary flow intensity and pattern at the inflow section: since the flow is fully developed in this section, i.e. it is obtained in the assumption of ∂w ∂z = 0, u and v form closed streamlines.

3. Analysis of the velocity profile at the inflow section to put in evidence the regions of high and low shear rate.

4. Temporal and spatial development of the degree of encapsulation.

From the comparison of the points 2 and 4, it should be possible to determine how predictive is the secondary flow pattern of a quasi2D simulation with respect to the actual three-dimensional encapsulation development. This point is fundamental since if both are related to each other, as it is expected to be, it means that, for a given pair of materials, the quasi2D analysis yields results predictive of the threedimensional flow but with a considerable reduction of the computational effort.

Newtonian fluid -Giesekus fluid pairing

The Newtonian/Giesekus combination offers the possibility to arbitrarily set the viscosity ratio while keeping the jump of N 2 across the interface always of the same available to the knowledge of the author and the fact that the radius of curvature of the corner blob is small when compared to the interface width and mesh size indicates that a further mesh refinement should be made in order for this result to be quantitatively relevant (see section (4.5)). However the latter is definitely consistent with the experimental observations reported by Mauffrey in [START_REF] Mauffrey | Étude numérique et expérimentale du phénomène d'enrobage dans les écoulements de coextrusion[END_REF] and summarized in section (1.1.3).

With respect to the steady state configuration, the encapsulation degree can be determined as:

E d = L 1 + L 2 + L 3 D (5.3)
Where L 1 , L 2 and L 3 are defined as in fig.

(5.12(b)). This definition is more appropriate than the classical definition relying on the fraction of perimeter wet by the less viscous fluid, because it includes the interface deviation from the midplane L 1 . This implies that even after the two contact lines have come to join, E d still increases as the encapsulated component moves towards the section center even if the section perimeter is completely wet by one fluid. As reported in [START_REF] Teixeira | Étude expérimentale et numérique du phénomène d'enrobage en coextrusion[END_REF][START_REF] Mauffrey | Étude numérique et expérimentale du phénomène d'enrobage dans les écoulements de coextrusion[END_REF], when plotted against z, as in fig.

(5.12(a)), the encapsulation rate can be approximately represented by a piecewise linear function with a discontinuous steepness at the point where the contact line meets the section corner.

The fact that, independently of the viscosity ratio, the encapsulation always occurs in the same sense and that the Newtonian fluid encapsulates the Giesekus fluid gives credit to the hypothesis that the phenomenon is intrinsically related to the sign of the jump of N 2 across the interface.

To further reinforce this result, all the simulations are repeated for the second case

Giesekus fluid -Giesekus fluid pairing

The Giesekus/Giesekus combination offers the possibility to observe the inversion of the encapsulation sense at high shear rates as first observed by Southern and Ballman in [START_REF]Additional observations on stratified bicomponent flow of polymer melts in a tube[END_REF][START_REF] Southern | Stratified bicomponent flow of polymer melts in a tube[END_REF]. This phenomenon was originally attributed to the presence of a viscosity crossover. In this section the inversion of encapsulation can be explained only as an effect of N 2 with arguments consistent with the ones provided in the previous section. For this purpose, the parameters of the Giesekus model are choosen as The values of the mobility parameter a, although a little too high when compared to real polymer melts, for which the ratio of the second to the first normal stress difference usually ranges between 0.1 and 0.3, are assumed in such a way as to magnify the recirculation effects and be able to observe a significant encapsulation within a relatively short die length.

λ 2 /λ 1 = 0.2, η 0 2 /η 0 1 = 0.
The same simulation parameters and geometry as in the previous section are retained From the analysis of the secondary flow pattern in the inflow section the following predictions can be made: In the first case, the flow of a less elastic and viscous fluid encapsulating a more elastic and viscous one yields analogous results as for the case of the Newtonian/Giesekus pairing. The fluid 2 is pushing down at the wall and it is being pushed up in the midplane.

In the second case, the less elastic fluid 2 is encapsulated by the more elastic fluid 1.

By inspecting the three-dimensional interface evolution in fig.

(5.17(c)) to (5.17(k)), it is possible to notice that the deformation in the midplane is not in the same sense as the one suggested by the analysis of the secondary flow direction in the inflow section. While the latter indicates a downwards displacement in the midplane section, the three-dimensional interface becomes slightly bulged upwards. This represents the only case where the results of the quasi2D flow with a fixed interface and the actual three-dimensional calculation differ from each other.

The third and most interesting case is the one at W i = 1.1. For this configuration the fluid 1 has a lower absolute N 2 than fluid 2 at the wall, where the shear rate is maximum, while the situation is reversed when moving along the interface towards the low shear region, i.e. the center of the section. The inversion of the N 2 jump sign along the interface produces a deformation that is characterized by a double curvature. The evolution of the interface, shown in fig. (5.18(c)) to (5.18(k)), is consistent with the observations on the shape interface between PS and HDPE reported by Han in [START_REF] Han | A study of bicomponent coextrusion of molten polymers[END_REF]. This particular shape, is sometimes referred to as "gullwing" interface, and it was usually attributed to the viscosity crossover.

If the same numerical experiments are repeated with two viscoelastic fluids presenting a viscosity crossover but not a second normal stress difference crossover, the inversion of the encapsulation sense is not observed. At this point the results seem to clearly point out that phenomenon of encapsulation is driven by the normal stresses in the cross section plane and not by a viscosity difference. However, for real fluids these two elements are not always separable and the fluids having higher viscosities in general are also characterized by higher normal stress differences. In the case of two-phase flow under the assumption of a static flat interface and zero derivatives in the flow direction, which is the case for the inflow section, the interface is not at equilibrium, as confirmed by the presence of a non zero normal velocity on it. The non equilibrium is given by the unbalanced normal stresses on both sides of the interface. The coalescence of the contour lines of τ xx -τ yy near the contact line, observable in fig. (5.14(b)) and (5.14(b)) for instance, indicates that in the limit of sharp interface there is a jump of the normal stress difference across the interface. If the analysis is restricted to the interface portion close to the wall, depicted in the schematics of fig. (5.19), where this stress difference corresponds to N 2 , following the same idea as in [START_REF] Yue | Elastic encapsulation in bicomponent stratified flow of viscoelastic fluids[END_REF], it is possible to conclude that the interface side facing the region of fluid having a greater negative N 2 , experiences a greater tensile stress in the second direction (here y) than the side facing the region of lower N 2 . This means that the net force exerted by the fluids on the interface tends to pull it upwards in the region occupied by the fluid possessing the greater absolute value of N 2 . This scenario is consistent with the direction of the flow near the contact line In three-layers systems the interface is shifted from its initial position h i = 0.25 even if the viscosity ratio and the flow rate ratio are 1. To obtain the value of the interface displacement, the mass conservation equation is applied to the control volume ABCD. The flow rates through AB, AD and BC are given by: Here the flow rate through the inflow section has been set to one fourth of the total flow rate Q, while the flow rates through the solid wall and the interface are zero since V • n is uniformly zero at these boundaries. To evaluate the flow rate through DC, the channel section at z = ∞, the velocity profile is assumed as fully developed.

Hence: effect. In the next section this phenomenon is analyzed in combination with the effects of viscosity for the first regime in three-dimensional flow.

Analysis of the first regime of encapsulation

The three-dimensional flow analyzed in this section is characterized by the same parameters as in section (5.2.1). The same simulation parameters are also retained with the only exception that for this case the flow has two planes of symmetry, so only one quarter of the domain is considered. Instead of separating the two cases A -B -A and B -A -B, they are treated as if the ratio of the core to the skin layer viscosities varied from 0.1 to 10. For instance, the cases η 1 /η 2 = 0.1 and η 1 /η 2 = 10 correspond to the same viscosity ratio but the former refers to the configuration A -B -A while the latter to B -A -B.

From the observation of the steady state interface shape in the symmetry plane of the channel in fig. (5.21), the it can be concluded that:

1. As for the two-dimensional case, for a viscosity ratio of 1 the core layer is squeezed by the two skin layers due to a pure geometric effect 2. This phenomenon is favored by increasing the viscosity of the skin layers.

3. The same initial layer repartition can be maintained by increasing the viscosity of the core layer. In the numerical simulations this seems to occur for a viscosity ratio between 5 and 10.

The evolution of the interface deformation in time is shown by fig. (5.22), for the configurations corresponding to η 1 /η 2 = 0.1 and η 1 /η 2 = 10. In the former case the initial perturbation grows in the flow direction forming two symmetric crests that fold up and merge again with core layer, while in the latter the effects of flow rates and viscosity compensate each other and thus the initial perturbation is damped indicating a greater stability of the flow in case the core layer is made of the more viscous fluid.

Analysis of the second regime of encapsulation

Since multi-layer coextrusion is used mainly for the production of thin films, the analysis of the second regime for the three-layers configuration is extended to rectangular sections too. The combination of materials considered is the same as the one for the case N EW T GIE1, with a viscosity crossover at λ H γ = 1. This choice is motivated by results of the section (5. effect of the quasi2D assumption as discussed in section (5.2.3), however the threedimensional evolution of the interface, shown in fig. (5.28), confirms that the core layer made of Newtonian fluid actually pushes into the skin layers. This effect is even more evident for W i = 0.5995 for which the overall encapsulation degree is smaller.

This result apparently undermines the validity of the N 2 hypothesis, however it should be noted that as the aspect ratio of the section is increased, the flow in the symmetry plane of the channel tends to the two-dimensional flow between two parallel plates. In this limiting case the second normal stress difference effects are canceled out and encapsulation also vanishes. Hence, for high aspect ratios the mechanism of the interface deformation in the symmetry plane is somewhat similar to the first regime mechanism. It has already been emphasized in section (5.1) how the first regime cannot be completely excluded from the second regime simulations due to the shear thinning behavior of the Giesekus fluids. In fact, as the encapsulation evolves in space the mean shear rate of two adjacent layer changes and consequently their viscosity ratio changes as well. This means that, since the flow rate ratio is fixed, the interface will be subject to a shift due to the change in the viscosity ratio and the higher the aspect ratio, the more this effect becomes predominant.

When the aspect ratio is high enough, the effect of geometry can so prevail that the core layer can break up halfway between the wall and the plane of symmetry and the Newtonian fluid remains entrapped in the middle leading to the anomalous encapsulation. This tendency is favored by incompressibility in quasi2D simulation, as seen in section (4.4), where the phenomenon occurs for a section aspect ratio of 4 : 1 while in the experiments of Khomami it is observed for aspect ratios greater than 10 : 1. dimensionless time units. This is because the point where they are supposed to meet lies in a symmetry plane where the normal velocity is zero. Hence, as the contact line approaches this region, it progressively slows down. The coalescence is due to short range interactions occurring when the two interfaces profiles overlap but this phenomenon is observable on an even longer time scale (due to the high P e) than that of secondary flows. Indeed, most of experimental results indicate a time of 15 to 30 minutes necessary to ensure the achievement of steady state conditions, which is a lapse of time one order of magnitude greater than the one of the numerical simulations presented in this chapter. Un facteur supplémentaire est représenté par la géométrie de la section: pour des rapports hauteur/largeur très petits le mouvement de l'interface dans le plan de symétrie de la section est dirigé dans le sens opposé à celui imposé par la différence de N 2 et, au delà d'une valeur critique de ce paramètre, une inversion du sens de l'enrobage peut être observée.

6 Conclusions and recommendations

General conclusions

Phase-field modeling is found to be an attractive approach to simulate bicomponent coextrusion flows as an alternative to classic interface-fitted methods. Although the developed code can only be considered as first version deserving improvements in order to achieve the computational performance required to obtain significant quantitative results, the overall method proved to be able to take into account for all the different mechanisms involved in the encapsulation phenomenon.

Using numerical schemes on uniform meshes that are only second order in space can be a limiting factor to the global accuracy of the model, since in order to correctly model the surface tension forces it is necessary to accurately resolve the hyperbolic tangent profile of the interface. This requirement could only be met for fairly wide interfaces or, alternatively, by using a higher order numerical scheme. In case of coextrusion flows, surface tension plays a minor role and satisfying results can be found even for thinner interfaces since the Capillary number is high enough to allow to neglect the interface tension.

The numerical solvers adopted in the present work to invert the discrete Stokes problem and the phase-field discrete equation are shown to perform adequately enough for the problem at hand, although the implementation of the Vanka smoothing and Non-linear Gauss Seidel techniques takes advantage of the relative simplicity of the geometry considered. The extension to more general unstructured meshes would be non trivial and it would require a considerable programming effort to try to make these solvers independent on the underlying geometry.

Even if the adopted schemes have an overall second order accuracy, only first order accuracy can be obtained when capturing a discontinuity in the flow field. The latter result is independent on the choice of the particular phase-field model and numerical scheme adopted since it is found from the examination of the particular case in which the phase field variable φ is just a passive scalar transported by the velocity field and using high order reconstruction for the treatment of the advection term. This accuracy limit justifies the presence of a diffuse interface introducing a first order error o (ǫ), since the price for controlling the interface behavior, and especially its diffusivity, is the same as the one due to the presence of a discontinuity in the numerical solution. On the other hand, this loss of accuracy near the interface might undermine the efficacy of the method itself.

In fact, from the analysis of the numerical solution of three-dimensional coextrusion flows, the second normal stress difference jump acting on the interface is found to be the real "engine" of encapsulation. Logically, this jump must be accurately resolved in order to get correct results. Unfortunately the encapsulation principal driving factor results from the resolution of the flow region where the accuracy of the method is the lowest, that is the interface. Another factor determining the accuracy of the solution in this region, is the adoption of the linear mixing rule; even though the application of more sophisticated laws, like reptation mixing rules (see [START_REF] Leonardi | Rheological models based on the double reptation mixing rule: The effects of a polydisperse environment[END_REF]) for the stress tensor or harmonic mixing rule for the density and the viscosity (see [START_REF] Liu | A phase field model for the mixture of two incompressible fluids and its approximation by a fourier-spectral method[END_REF][START_REF] Khatavkar | Diffuse interface modeling of droplet impact[END_REF]), is not guaranteed to produce better results since the phase diffusion in in the phase-field method is just a numerical artifact used to control the interface behavior while the real mixing occurs at much smaller scales.

With in mind the limits exposed above, the phase-field model is still useful for the comprehension of the physics behind the encapsulation phenomenon. The inspection of the flow solutions obtained under the assumption of quasi2d approximation and static (fixed) interface allowed an estimate of the Deborah number of encapsulation. Results show that, independent of the Weissemberg number of the flow, the Deborah number is at most of the order of magnitude of 10 -3 for square ducts. This value represents an overestimate since it makes use of the initial value of the average cross sectional velocity to estimate the characteristic encapsulation time. As encapsulation advances in time, the flow evolves towards an axisymmetric configuration (full encapsulation) where the secondary flow intensity tends to zero. This indicates that the Deborah number is almost zero and that elasticity has little to do with the layer rearrangement observed, since coextrusion is an almost viscosimetric flow. So, instead of discriminating viscous encapsulation from elastic encapsulation as usually presented by the specialized literature, these are interpreted as two manifestations of the fluid viscoelasticity. In the intent of characterizing this phenomenon, the first and the second regime of encapsulation are examined independently. The first regime, consisting of a rapid interface shift, depends on the pressure gradient conditions with whom the two polymeric fluids come into contact. The layer possessing the greater pressure gradient pushes into the one presenting the lowest one until the equilibrium is restored. This phenomenon occurs with little interface deformation and depends on the layers repartition, flow rate and viscosity. In two-layers coextrusion, for equal flow rates and layer thicknesses, the first regime translates into the protrusion of more viscous fluid into the less viscous one. In three-layers configurations of the type A -B -A for balanced layer repartition and flow rates (i.e 0.25/0.5/0.25), the first regime is complicated by the fact that interfaces do not lie in the section symmetry plane and even if the viscosity ratio is 1 the pressure gradient imbalance provokes the squeezing of the central layer. For more general configurations both these two effects are observed.

The second regime is due to the unbalanced normal stress difference acting on the surface. The interface reacts to this unbalance by moving towards the fluid having the greatest absolute N 2 near the wall and towards the one having the smallest N 2 in vicinity of the midplane, thus determining an interface shape that is convex on the side of the fluid with the greatest negative N 2 . If the two components have a crossing point in the curve of Ψ 2 then, under particular shear rate conditions, the interface assumes the so called "gull-wing" shape, consisting of a double curvature.

The flow direction, and consequently the interface deformation, is not always consistent with the direction suggested locally by the N 2 criterion discussed above, because it often results from the competition of the first and second regimes combined together. These two are non-linearly related and thus they are not simply superposable. In fact, as the interface deforms for effect of the second regime, the viscosity in each component varies due to the readjustment of the velocity field to the new conditions and the resulting change in the shear rate. A crucial role is played by the geometry. For high values of the aspect ratio of the rectangular section, the first regime effects prevail in center of the section while the second regime effects become more evident at the wall. In three-layers coextrusion, for the case of a central layer having a smaller N 2 than the skin layers, the interface moves towards the upper and lower walls in the symmetry plane of the section, opposed to what it is dictated by the N 2 jump sign. Consequently, if the central layer is sufficiently thin it will rupture on both sides and become surrounded of the fluid having the greatest negative N 2 . This phenomenon is referred to as anomalous encapsulation.

Recommendations

The main difficulties encountered in the simulation coextrusion flows by the phasefield method are due to the loss of accuracy of the model in the region near the interface. Two solutions can be envisaged to fix this problem. The first consists in modifying the phase field equation in order to include additional terms that are specifically tailored to eliminate or limit some known issues of the model. Since the chemical potential terms can also be related to the curvature of a diffuse interface having a hyperbolic tangent profile, the Cahn-Hilliard equation can be interpreted as curvature-driven flow equation in which a free surface evolves as to minimize its global curvature. This means that the Cahn-Hilliard model is not shape-preserving unless a very high resolution is used (i.e. the interface width to radius of curvature ratio must tend to zero). In practice, since the interface is artificially diffused, even in absence of an external flow and surface tension effects, a closed surface of any shape evolves into a sphere. For this reason, an additional term is included in order to cancel the curvature leading order effect in the equation in such a way that the phase-field diffusion only acts in the direction normal to the interface while the shape of the latter is entirely determined by the velocity field and capillary effects are only due to the surface tension term in the Navier-Stokes equation. The resulting phasefield equation is no longer conservative and that is the reason why this solution is often applied to the Allen-Cahn equation since it is already non conservative but only second order instead of fourth order. This method was first introduced by Floch in [START_REF] Folch | Phase-field model for Hele-Shaw flows with arbitrary viscosity contrast. I. Theoretical approach[END_REF][START_REF]Phase-field model for Hele-Shaw flows with arbitrary viscosity contrast. II. Numerical study[END_REF], while a similar approach, called the Advected-Field method can be found in [START_REF] Biben | An advected-field approach to the dynamics of fluid interfaces[END_REF]. Recent developments lead to hybrid level-set/phase-field models relying on a modified Cahn-Hilliard equation and a sharp interface reconstruction like in [START_REF] Sun | Sharp interface tracking using the phase-field equation[END_REF].

The second method is to use Adaptive Mesh Refinement (AMR) to improve the resolution only in the neighborhood of the interface [START_REF] Martin | A cell-centered adaptive projection method for the incompressible Navier-Stokes equations in three dimensions[END_REF]. The attractiveness of AMR is in that the Cahn-Hilliard equation could be retained in its original thermodynamic sense without any loss of accuracy and the simulation could take into account for both the macroscopic flow and the phase intermixing occuring at the interface scales.

Although this alternative would be closer to the real physics of diffuse interfaces, the macroscopic and the microscopic scales are so widely separated that the direct simulation of both is still an out-of-reach goal. Nonetheless, in most practical cases the phase-field does not need a full scale simulation to converge, so AMR still represents a completely viable option. Another promising feature of the phase-field method is that it can be efficiently applied to multi-phase systems in a straightforward way [START_REF] Kim | Phase field computations for ternary fluid flows[END_REF][START_REF]Phase-Field Models for Multi-Component Fluid Flows[END_REF].

Once the two-phase flow model improvement is assessed, the numerical study of encapsulation could be extended some other aspects including:

1. the influence of the wall affinity and the choice of a more suitable wall condition for the contact line; 2. the effect of the surface tension, which should be hindering the second regime; 3. the effect of a variable section die, for instance a tapered die; 4. the effect of a multi-mode constitutive equation. The latter point is particularly important since the point of maximum in the curve of the secondary flow intensity seems to be associated with the presence of a single relaxation time. There is only one Weissemberg number condition for which the recirculation and thus the coextrusion are maximized. Since real polymer melts posses a spectrum of relaxation times, a more realistic multi-mode model is expected to predict a larger interval of operational conditions in which the manifestation of encapsulation is more substantial.

Conclusions et perspectives

1. 1 4 1. 1 . 2 A 1 vof Figures 1 . 1

 1412111 Review of the experimental studies . . . . . . . . . . . . . . . . . . . 4 1.1.1 1972-1976: The circular die . . . . . . . . . . . . . . . . . . . dual mechanism hypothesis . . . . . . . . . . . . . . . . . . 1.1.3 1990-2002: The rectangular die . . . . . . . . . . . . . . . . . 1.2 Theoretical Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 The minimum viscous dissipation principle . . . . . . . . . . . 1.2.2 The first regime of encapsulation . . . . . . . . . . . . . . . . 1.2.3 The second regime of encapsulation . . . . . . . . . . . . . . . 1.3 Numerical Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Aim of the present work . . . . . . . . . . . . . . . . . . . . . . . . . The phase-field model . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Derivation of the governing equations . . . . . . . . . . . . . . 2.1.2 Physics of the phase-field model . . . . . . . . . . . . . . . . . 2.1.3 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . List Schematic of the layer rearrangement occurring in two-layers coextrusion in a circular die. . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Shear viscosity of polystyrenes A nd B (a). Interface profiles for different shear rates (b). Southern and Ballman, 1973 [49]. . . . . . . 1.3 Shear viscosity of LDPE and PS (a). Schematic of the experimental appartus (b). White et al. 1973 [26] . . . . . . . . . . . . . . . . . . 1.4 Material functions of HDPE, LDPE, PS and PMMA. Shear viscosity (a). First normal stress difference (b). Second normal stress difference (c). Lee and White, 1974 [36] . . . . . . . . . . . . . . . . . . . . . . 1.5 Interface shape for HDPE/LDPE extrudates as a function of the capillary die length and time. L/D = 20 (a). L/D = 10 (b). L/D = 5 (c). [36]. Lee and White, 1974 [36] . . . . . . . . . . . . . . . . . . . 1.6 Percentage of the wall perimeter wet by the most viscous fluid as a function of the viscosity ratio (a) and capillary length (b). Lee and White, 1974 [36] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Shear viscosity for HDPE, LDPE and PS. Han, 1973 [22]. . . . . . . 1.8 Interface profiles for the PS-LDPE configuration for L/D = 4. Q P S = 10.9cc/min and Q LDP E = 9.9cc/min (a). Q P S = 8.3cc/min and Q LDP E = 53.8cc/min (b). Han, 1973 [22]. . . . . . . . . . . . . . . . 1.9 Interface profiles for the PS-HDPE configuration for L/D = 18. Q P S = 8.3cc/min and Q HDP E = 11.3cc/min (a). Q P S = 8.3cc/min and Q HDP E = 54.2cc/min (b). Han, 1973 [22]. . . . . . . . . . . . . 1.10 Extrudate bending towards the high viscosity component occurring at the exit of the die. Everage, 1975 [1] . . . . . . . . . . . . . . . . . Bending angle and interface profile as a function of the capillary length L/D. Everage, 1975 [1] . . . . . . . . . . . . . . . . . . . . . . 1.12 Interface deformation for a two-layer polystyrene structure as it flows down a square section die. The initial repartition is 20 : 80. Dooley, 2002 [12] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.13 Coextrusion of two layer of PS 7240 at different temperatures. ∆T = 15 • C (a). ∆T = 30 • C (b). Mauffrey, 1998 [39] . . . . . . . . . . . . . 1.14 Coextrusion degree (percentage of wet perimeter) of the pairing PS 1160 and PS 7240 for different aspect ratios of the rectangular section and surface roughness. Aspect ratio = 1 (a). Aspect ratio = 2 (b).

1. 15 2 . 7

 1527 Coextrusion degree for the polyethylene-polyamide pairing evidencing the presence of a plateau.Mauffrey, 1998 [39] . . . . . . . . . . . . .1.16 Stratified flow between parallel plates. Side-by-side configuration (a). Sheath-core configuration (b). . . . . . . . . . . . . . . . . . . . . . . 1.17 Pressure gradient as a function of reciprocal viscosity ratio for the side-by-side and sheath-core configurations. Flows are considered at equal flow rates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.18 Interface profile as a function of the channel length, obtained by solving the differential equation derived by Everage [1] for different viscosity ratios. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.19 Comparison between an interface-fitted mesh, Sizaire, 1998 [48] (a) and a diffuse interface mesh, Yue, 2006 [65] (b). . . . . . . . . . . . 2.1 Evolution in time of spinodal decomposition. Simulation is carried out on a spatially periodic domain box. . . . . . . . . . . . . . . . . . 2.2 Bulk free energy density as a function of the phase composition. . . . 2.3 One dimensional equilibrium solution of the phase-field model. . . . . 2.4 Slipping mechanism of a diffuse interface contact line. . . . . . . . . . 2.5 Shear viscosity as a function of shear rate for shear-thinning fluids. . 2.6 Shear viscosity of single mode Giesekus model for different values of the non linear parameter a and η s = 0. . . . . . . . . . . . . . . . . . List of Figures xiii First and second normal stress difference coefficients of single mode Giesekus model for different values of the non linear parameter a and η s = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Numerical molecule for the Marker and Cell method. V and S denote the locations of the velocity components and strain rate tensor components. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Definition of the stencil for the TVD interpolation of φ. . . . . . . . . 3.3 Velocity red-black decoupling occurring for staggered arrangement of velocity and stress fields. . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Convergence history of the non linear GS/NGMRES method applied to the phase-field problem. . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Schematics of the restriction and prolongation operations on a staggered grid. Cell-centered residual in the node I is obtained by averaging fine grid values of nodes A to H (left). Staggered residual component in the node E is obtained by averaging the fine grid values corresponding to nodes A to D . . . . . . . . . . . . . . . . . . . 3.6 Schematic of a subdomain block. face, edge and corner ghost cells are represented only on one side for simplicity. . . . . . . . . . . . . 3.7 Speed-up factor for the parallel solution of the SINKER problem described in section ((4.2)) for a 256 3 grid (non balanced) and a 384 3 grid (balanced). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Comparison between analytical and numerical solution for Poiseuille flow of a Newtonian fluid in a square section. Velocity profiles are plotted at equally distanced positions along y (a). Numerical error as function of the mesh resolution showing the second order accuracy of the scheme (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . .4.2 Comparison between analytical and numerical solution for the Oldroyd-B flow between parallel plates with β = 0.1. (a) Velocity profile; (b) Shear stress; (c) Normal stress. (d) phase field. . . . . . . . . . . . . . 4.3 Error convergence for the normal stress and the phase field. . . . . . . xiv List of Figures 4.4 Schematic of the SINKER problem also showing streamlines for ∆η = 10 -3 . The central block in red is made of the more viscous and more dense fluid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Convergence history of the SINKER problem solved with GCR+FVS on 64 × 64 × 64 grid with a viscosity jump defined over 3.5 cells and different values of the viscosity contrast. . . . . . . . . . . . . . . . . 4.6 Solver sensitivity to h, for ∆η = 10 -3 and a viscosity jump defined over 3.5 cells. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Schematic of the domain used for the calculation of secondary flow in square ducts. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Intensity of secondary flows in a square duct scaled by the intensity of the main flow (a), and by L/λ H (b). Comparison with numerical data of Yue et. al. [66, fig. 3] . . . . . . . . . . . . . . . . . . . . . . 4.9 Stress difference τ eyy -τ exx for De = 1, a = 0.5 β = 0.1. (a) Calculation by Yue et al. [66, fig. 2(c)] (b) Present simulation. (c) Streamlines of secondary flow. . . . . . . . . . . . . . . . . . . . . . 4.10 Schematic of the domain used for quasi-2d calculation. The aspect ratio of the rectangular section is equal to 4. . . . . . . . . . . . . . . 4.11 Temporal evolution of anomalous encapsulation (rescaled time). (a) t = 0 . (b) t = 27.85 . (c) t = 50.13 . (d) t = 83.75 . (e) t = 108.61 . (f) t = 136.46 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.12 Comparison between the numerical solution in fig [65, fig. 10(a)] at time t = 76.21 (bottom) and the present solution at t = 136.46 (top). Temporal evolution of the interface is slower due to the secondary motions being weaker because of the greater Deborah number in our simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 Shear viscosity material function of the Giesekus model for PS 1161 and PS 4801 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 Schematic of the domain used for three-dimensional calculation of encapsulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15 Schematic of the three-layers coextrusion configuration. . . . . . . . . List of Figures xv 4.16 Mid-plane view of three layers coextrusion of Newtonian fluids with R η = 0.5 and Q 1 /Q 2 = 1 (flow rate ratio) at t = 0.25. (a) h = 0.03125; (b) h = 0.015625; (c) h = 0.0078125. . . . . . . . . . . . . . 4.17 Mid-plane view of three layers coextrusion of Newtonian fluids with R η = 0.1 and Q 1 /Q 2 = 5.25 (flow rate ratio) at t = 0.5. (a) h = 0.03125; (b) h = 0.015625; (c) h = 0.0078125. . . . . . . . . . . . . . 4.18 Temporal evolution of the PS 4801 layer. (a) t = 0, (b) t = 5, (c) t = 10, (d) t = 35, (e) t = 60. . . . . . . . . . . . . . . . . . . . . . . 4.19 Interface profiles for configurations with (right) and without (left) the separation plate. (a,e) z = 0.5, (b,f) z = 2.5, (c,g) z = 5, (d,h) z = 8 4.20 (a) Experimental interface profiles from z = 0 to z = 24. Data courtesy of Teixeira, 1996 [44] (b) Comparison between the numerical encapsulation at z = 5 and the experimental data at z = 5.54 . . . . 5.1 Schematic side view of the computational domain used for the simulations of the first regime. . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Comparison of two-layers coextrusion flows: η 2 /η 1 = 0.2 (left), η 2 /η 1 = 0.05 (right). t = 0 (a,e). t = 2 (b,f). t = 4 (c,g). t = 6 (d,h). . . . . 5.3 Interface profile and velocity profiles in the symmetry plane. η 2 /η 1 = 0.2 (left), η 2 /η 1 = 0.05 (right). t = 1 (a,e). t = 2 (b,f). t = 3 (c,g). t = 4 (d,h). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Interface profile in the cross section plane z = 1 at t = 6. η 2 /η 1 = 0.2 (a). η 2 /η 1 = 0.05 (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Interface profiles in the symmetry plane for different values of the viscosity ratio. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Interface position and entry length as functions of the viscosity ratio (a). Dimensionless pressure gradient and viscous dissipation as functions of the viscosity ratio (b). . . . . . . . . . . . . . . . . . . . 5.7 Evolution of the interface shape due to a sinusoidal flow rate perturbation in a two-layers stratified flow with Re = 10 -6 , η 2 /η 1 = 0.1. t = 4 (a). t = 5.5 (b). t = 7 (c). t = 8.5 (d). t = 10 (e). t = 11.5 (f) . xvi List of Figures 5.8 Lateral view of the computational domain used for the simulations of the second regime of encapsulation (a). Local coordinate system for the definition of N 2 at the wall (b) . . . . . . . . . . . . . . . . . . . 5.9 Shear viscosity as a function of the shear rate for Newtonian and Giesekus fluids. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 Intensity of the secondary flow as a function of the Weissemberg number for the configuration N EW T GIE1 (a). Visualization of the secondary flow in the left half of the square cross section. Newtonian fluid is red (upper) and Giesekus fluid is blue (lower). W i = 0.0278 (b). W i = 0.2154 (c). W i = 1.6681 (d). . . . . . . . . . . . . . . . . 5.11 Case N EW T GIE1 at W i = 0.2154. Inlet velocity profile at indicated x sections (a) and normal stress (b). Visualization of the deformation of the Newtonian fluid layer. t = 0 (c). t = 7.5 (d). t = 15 (e). t = 22.5 (f). t = 30 (g). t = 37.5 (h). t = 45 (i). t = 52.5 (j). t = 60 (k). . . . . . . . . 5.12 Encapsulation level for the case N EW T GIE1 at W i = 0.2154 (a), as defined in the schematic (b). . . . . . . . . . . . . . . . . . . . . . 5.13 Intensity of the secondary flow as a function of the Weissemberg number for the configuration N EW T GIE2 (a). Visualization of the secondary flow in the left half of the square cross section. Newtonian fluid is red (upper) and Giesekus fluid is blue (lower). W i = 0.0278 (b). W i = 0.2154 (c). W i = 1.6681 (d). . . . . . . . . . . . . . . . . 5.14 Case N EW T GIE2 at W i = 0.2154. Inlet velocity profile at indicated x sections (a) and normal stress (b). Visualization of the deformation of the Newtonian fluid layer. t = 0 (c). t = 5.25 (d). t = 10.5 (e). t = 15.75 (f). t = 21 (g). t = 26.25 (h). t = 31.5 (i). t = 37.75 (k) . . . . . . . . . . . 5.15 Case GIEGIE2. Shear viscosity (a). Second normal stress difference coefficient (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.16 Intensity of the secondary flow as a function of the Weissemberg number for the configuration GIEGIE2 (a). Visualization of the secondary flow in the left half of the square cross section. Fluid 1 is in red (upper) and fluid 2 fluid is in blue (lower). W i = 0.2154 (b).W i = 1.1 (c). W i = 3 (d). . . . . . . . . . . . . . . . . . . . . . . . . List of Figures xvii 5.17Case GIEGIE2 at W i = 3. Inlet velocity profile at indicated x sections (a) and normal stress (b). Visualization of the deformation of the fluid 2 layer. t = 0 (c). t = 7.5 (d). t = 15 (e). t = 22.5 (f). t = 30 (g). t = 37.5 (h). t = 45 (i). t = 52.5 (j). t = 60 (k). . . . . . . . . . . . . . . . . . . 141 5.18 Case GIEGIE2 at W i = 1.1. Inlet velocity profile at indicated x sections (a) and normal stress (b). Visualization of the deformation of the fluid 2 layer. t = 0 (c). t = 7.5 (d). t = 15 (e). t = 22.5 (f). t = 30 (g). t = 37.5 (h). t = 45 (i). t = 52.5 (j). t = 60 (k). . . . . . . . . . . . . . . . . . . 142 5.19 Schematic of mechanism driving the interface deformation at the wall for the cases N EW T GIE1 and N EW T GIE2. . . . . . . . . . . . . . 143 5.20 Schematic of the three-layers configuration. . . . . . . . . . . . . . . . 145 5.21 Three-layers configuration: comparison of the interface position in the symmetry plane for different viscosity ratios. . . . . . . . . . . . . 146 5.22 Comparison of three-layers coextrusion flows (only half channel): η 2 /η 1 = 0.1 (left), η 2 /η 1 = 10 (right). t = 0.75 (a,e). t = 1.5 (b,f). t = 2 (c,g). t = 3.08 (d,h). . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 5.23 Secondary flow intensity as a function of the Weissember number for rectangular and square sections with A -B -A and B -A -B configurations. Non-dimensionalization with respect to the main flow intensity W m (a), and λ H /D (b). . . . . . . . . . . . . . . . . . . . . 149 5.24 Secondary flow for the three-layers configuration in the square section at W i = 0.2154 (maximum recirculation intensity). A -B -A (a). B -A -B (b). B is the Newtonian fluid and A is the Giesekus fluid. 150 5.25 Secondary flow for the three-layers configuration in the rectangular section at the Weissember numbers yielding the maximum recirculation intensity. A -B -A and W i = 0.2154 (a). B -A -B and W i = 0.5995 (b). B is the Newtonian fluid and A is the Giesekus fluid.151 5.26 Interface deformation for the three-layers configuration A -B -A with square section at W i = 0.2154. t = 7.5 (a). t = 15 (b). t = 22.5 (c). t = 1 IntroductionMultiphase confined flow is one of the most investigated research topics because of its numerous industrial applications. The most interesting case is without doubt the flow of stratified oil and water layers. The relevance of this system is given by the tendency of the water to migrate toward the pipe wall in reason of its lower viscosity, giving rise to the so called core-annular flow (CAF). This particular configuration is characterized by a lower pressure drop than the single component flow produced by the same oil flow rate. The simple observation that the addition of a 8 -10% of water in crude oil pipelines can considerably reduce the pressure gradient needed to sustain the flow, leaded to the registration of several patents at the beginning of the 20 th century. Consequently, a significant research effort has been focused on the study of the stratified bicomponent flows, the proliferation of different regimes (slug flow, plug flow, CAF...) and their domain of stability. CAF is extensively documented in the remarkable work of Joseph who defined it as a gift of nature,[START_REF] Joseph | Core-annular flows[END_REF].Another technological application of paramount interest is represented by coextrusion. This process allows to combine several materials, in most cases polymers, into a single product in order to obtain specifically tailored properties otherwise impossible to attain with a single component. A noteworthy example of a coextruded product is represented by conjugate fibers. These are produced by coextruding two components in a side-by-side semicircular configuration. The difference of the thermal expansion coefficients of the materials causes the buckling of the fiber upon solidification at the exit of the spinneret, giving it a strong resemblance to natural wool fibers.

Figure 1 . 1 :

 11 Figure 1.1: Schematic of the layer rearrangement occurring in two-layers coextrusion in a circular die.

  Apart from flow instabilities, which represent a restrictive factor to the operational conditions range, polymer coextrusion flows are interested by unique phenomena that are due to the viscoelastic nature of the polymer melts. More specifically, when two different polymers are coextruded in a common die, the originally stratified flow can evolve toward a configuration where one component completely encapsulates the other in reason of their different rheological properties as shown in fig. (1.1). Although viscoelastic encapsulation is similar to the CAF transition in the sense that the less viscous component generally encapsulates the more viscous one, its driving mechanisms are basically different and thus the same considerations that are uniformly accepted for Newtonian fluids cannot be straightforwardly applied to this phenomenon. Experimentation offers a solid approach to investigate the subject of the interface rearrangement occurring in stratified flows of viscoelastic fluids; however the interpretation of the experimental results must rely on theoretical basis which are not as well established as for Newtonian fluids. Numerical simulation can offer a useful tool to complement and support the experimental observations.

  1972 and 1976 a number of papers were published on the subject of coextrusion in circular dies because of the increasing interest raised by the advent of this kind of technology in the textile industry. The interest on the subject was renewed again at the distance of twenty years and studies on rectangular dies were carried on in the attempt of understanding the physics of stratified flows characterizing the production of thin multilayer films and possibly finding a cure to the irregularities in the interface configuration observed in such flows. It seems only logical to review these periods separately, especially considering the technological gap between both and the advances in the understanding of the rheology of polymer melts achieved within two decades. Without any claim of exhaustiveness on the matter, the principal results and observations are presented in two of the three sections that follow while a third is dedicated specifically to the experiments that led to the identification of two different encapsulation regimes.

Figure 1 . 2 :Figure 1 . 3 :Figure 1 . 4 :

 121314 Figure 1.2: Shear viscosity of polystyrenes A nd B (a). Interface profiles for different shear rates (b). Southern and Ballman, 1973 [49].

Figure 1 . 5 :Figure 1 . 6 :

 1516 Figure 1.5: Interface shape for HDPE/LDPE extrudates as a function of the capillary die length and time. L/D = 20 (a). L/D = 10 (b). L/D = 5 (c).[36]. Lee and White, 1974[START_REF] Biing | An Experimental Study of Rheological Properties of Polymer Melts in Laminar Shear Flow and of Interface Deformation and Its Mechanisms in Two-Phase Stratified Flow[END_REF]

Figure 1 . 7 :

 17 Figure 1.7: Shear viscosity for HDPE, LDPE and PS. Han, 1973 [22].

Figure 1 . 8 :

 18 Figure 1.8: Interface profiles for the PS-LDPE configuration for L/D = 4. Q P S = 10.9cc/min and Q LDP E = 9.9cc/min (a). Q P S = 8.3cc/min and Q LDP E = 53.8cc/min (b). Han, 1973 [22].

Figure 1 . 9 :

 19 Figure 1.9: Interface profiles for the PS-HDPE configuration for L/D = 18. Q P S = 8.3cc/min and Q HDP E = 11.3cc/min (a). Q P S = 8.3cc/min and Q HDP E = 54.2cc/min (b). Han, 1973 [22].

Figure 1 . 10 :

 110 Figure 1.10: Extrudate bending towards the high viscosity component occurring at the exit of the die. Everage, 1975 [1]

Figure 1 . 11 :

 111 Figure 1.11: Bending angle and interface profile as a function of the capillary length L/D. Everage, 1975 [1]

Figure 1 . 12 :

 112 Figure 1.12: Interface deformation for a two-layer polystyrene structure as it flows down a square section die. The initial repartition is 20 : 80. Dooley, 2002 [12]

Figure 1 . 13 :

 113 Figure 1.13: Coextrusion of two layer of PS 7240 at different temperatures. ∆T = 15 • C (a). ∆T = 30 • C (b).Mauffrey, 1998 [39]

  (1.13(a)) and(1.13(b)). Near the wall, the equilibration of the temperature is more rapid. This means that the jump in the fluid properties across the interface is nil at the wall (i.e the coextruded fluids are the same) while it is non-zero in the inner region of the die. The fact that the flow still manifests a layer rearrangement suggests that the mechanism of the encapsulation is not a wall effect, but it is determined by the whole flow field. Encapsulation is only more evident near the boundaries because the main flow becomes progressively weaker approaching this region and so the residence time of the fluid particles flowing near the wall becomes of the same order of the encapsulation characteristic time. When the temperature difference is less than 5 • C no significant deformation in the interface is observed which justifies the hypothesis of isothermal flow within a tolerance of ±3 • C.For bicomponent coextrusion the frozen samples are extracted from within the die and sliced at equally spaced positions along the axes. In this way the evolution of coextrusion along the die is measured. The degree of encapsulation, defined as the perimeter fraction wet by the less viscous fluid, is well represented by a piecewise linear function of the distance from the die entrance (fig. (1.14(a)) to (1.14(d))). The discontinuity of the steepness occurs for the position at which the contact line meets the section corner. Mauffrey explained the presence of the discontinuity as a result of the extra energy needed by the contact line to turn around the corner. Indeed this energetic consideration should explain why for the case of polyethylene-polyamide pairing (PE-PA), the contact line stagnates in corner until the accumulated energy is enough to allow a further advancement (fig. (1.15)

Figure 1 . 14 :

 114 Figure 1.14: Coextrusion degree (percentage of wet perimeter) of the pairing PS 1160 and PS 7240 for different aspect ratios of the rectangular section and surface roughness. Aspect ratio = 1 (a). Aspect ratio = 2 (b).

Figure 1 . 16 :

 116 Figure 1.16: Stratified flow between parallel plates. Side-by-side configuration (a). Sheath-core configuration (b).

Figure 1 . 17 :

 117 Figure 1.17: Pressure gradient as a function of reciprocal viscosity ratio for the sideby-side and sheath-core configurations. Flows are considered at equal flow rates.

Figure 1 . 18 :

 118 Figure 1.18: Interface profile as a function of the channel length, obtained by solving the differential equation derived by Everage [1] for different viscosity ratios.

Figure 1 . 19 :

 119 Figure 1.19: Comparison between an interface-fitted mesh, Sizaire, 1998 [48] (a) and a diffuse interface mesh, Yue, 2006 [65] (b).

  revealed the coexistence of two different mechanisms: a first regime consisting of a rapid interface shift due to the different viscosities and flow conditions at which two polymer melt streams come in contact with each other, and a second regime consisting in a slow encapsulation of one layer around the other depending on the rheological properties of the components. Although in most of the experiments reported it is the less viscous fluid which encapsulates the more viscous one, a simple viscosity difference criterion fails to explain the mechanism of the second regime since it is not able to justify certain experimental observations like, for instance, the absence of encapsulation for Newtonian fluids. The theoretical analysis results point out that the second normal stress difference is behind the mechanism of the interface distortion.

34 Chapter 1 .

 341 IntroductionStokes itératif à viscosité variable capable de gérer des sauts de viscosité arbitraires entre les deux fluides et l'adoption d'un schéma semi-implicite pour l'intégration en temps des équations. Ce choix permet d'un côté de résoudre séparément les champs de vitesse, phase et contrainte, contrairement au schéma implicite adopté par Yue qui implique une solution couplée de ces équations, et en même temps de ne pas avoir de contraintes importantes sur la réduction du pas de temps pour la stabilité numérique comme pour la méthode explicite adoptée par Mompean et al.[START_REF] Mompean | Finite volume simulation of viscoelastic flows in general orthogonal coordinates[END_REF].This chapter illustrates the mathematical model for the treatment of coextrusion flows. The phase-field model for two-phase incompressible flow is derived and the surface tension forces modeling is discussed as well as the appropriate boundary conditions to be assigned. A family of differential constitutive equations is briefly presented. Finally equations are made non-dimensional and dimensionless parameters governing the flow are derived.

Figure 2 . 1 :

 21 Figure 2.1: Evolution in time of spinodal decomposition. Simulation is carried out on a spatially periodic domain box.

Figure 2 . 2 :

 22 Figure 2.2: Bulk free energy density as a function of the phase composition.

16 )

 16 As shown by figure(2.3), the function above expresses the fact that the phase composition varies across the interface with an hyperbolic tangent profile.The presence of the group β 2α in the hyperbolic tangent leads to the definition of natural length scale associated to the interface thickness ǫ.

Figure 2 . 3 :

 23 Figure 2.3: One dimensional equilibrium solution of the phase-field model.

  CH equation is fourth order in space and thus requires two boundary conditions due to the presence of the biharmonic operator. The assignment of proper boundary condition at the wall is non trivial and requires a certain understanding of the dynamics of the contact line motion. A first condition comes from the notion that the phase composition cannot vary through a solid wall due to diffusion. This yields a Neumann condition on the chemical potential.

Figure 2 . 4 :

 24 Figure 2.4: Slipping mechanism of a diffuse interface contact line.

  flow and Newtonian fluids, both equation(2.14) and condition(2.24) and can be formally derived from the minimization with respect to the fluid velocity of the functional associated to the sum of the kinetic energy diffusive dissipation and the rate of change of free energy[START_REF] Tiezheng Qian | A variational approach to moving contact line hydrodynamics[END_REF], thus supporting the interpretations of interface deformation phenomena as the result of the minimization of energy dissipation.It is easy to see that, for single component flows, the GNBC reduces to the standard Navier Boundary Condition (NBC). The difference between the two resides in the presence of an unbalanced Young stress term L (φ) ∂φ ∂t expressing the "elastic" force exerted by the interface on the flow. This term is the slip velocity equivalent of the distributed surface tension forces in the NS equations for the interior flow field.When the fluid is Non-Newtonian the contribution -η ∂vt ∂n + ∂vn ∂t in equation(2.24)

ForFigure 2 . 5 :

 25 Figure 2.5: Shear viscosity as a function of shear rate for shear-thinning fluids.

Figure 2 . 6 :

 26 Figure 2.6: Shear viscosity of single mode Giesekus model for different values of the non linear parameter a and η s = 0.

Figure 2 . 7 :

 27 Figure 2.7: First and second normal stress difference coefficients of single mode Giesekus model for different values of the non linear parameter a and η s = 0.

Equations ( 2 .

 2 14) together with (2.26) constitute a closed system supplied with boundary conditions discussed in section (2.1.3). The family of constitutive equations analyzed are hyperbolic, and Dirichlet boundary conditions are needed only on boundaries where V • n > 0.

  .41) is written for each mode and for both components. For formal simplicity the subscripts 1, 2, indicating the bulk phases, have been dropped. Deborah number De i is different from the global Deborah number De defined on the basis of the natural relaxation time of the most viscous component. Its definition varies for each mode and fluid component.

  de fluides viscoélastiques. L'hypothèse de base de la méthode de champ de phase est de pouvoir décrire un fluide binaire par une fonction φ représentant sa composition locale. L'interface entre deux fluides est donc remplacée par une région très mince où les deux phases sont "miscibles". Cette hypothèse de miscibilité locale se traduit par l'introduction d'un terme de diffusion non-linéaire dans l'équation hyperbolique de transport de phase. Ce terme est essentiel pour la régularisation de cette équation et le traitement de la discontinuité. En même temps il est dérivé de la théorie cinétique de séparation de phase proposée par Cahn and Hilliard[START_REF] Cahn | Free Energy of a Nonuniform System. I. Interfacial Free Energy[END_REF].Malgré le fait que l'interface soit résolue à une échelle typiquement plus grande que dans la réalité, ce type de modélisation ayant une base physique se prête naturellement à décrire les phénomènes de rupture, de coalescence d'interface et de glissement du point triple de contact entre l'interface et une paroi solide, ainsi qu'à neutraliser les effets de diffusion numérique dérivant de la discrétisation du terme convectif de l'équation. Après adimensionnement des équations, il apparaît que la physique d'une interface diffuse est régie par deux paramètres: le nombre de Cahn Cn, représentant la séparation de l'échelle microscopique où les effets de mélange entre les deux phases ont lieu, avec l'échelle macroscopique de l'écoulement, et le nombre de Péclet P e, indiquant le rapport entre les effets de transport convectif de phase et les effets de diffusion. Le choix de Cn doit respecter un critère de résolubilité du profil d'interface sur un maillage à pas uniforme fixé h qui impose une limite de 3.5h pour l'épaisseur de l'interface. De la même façon, le P e doit garantir un équilibre local entre le transport convectif et la diffusion de phase pour que l'interface puisse garder une épaisseur constante indépendamment de l'action compressive ou dilatante exercée par le champ de vitesse. Cela se traduit par la relation P e × Cn = 1.

  the discretized equations is based on standard programming language C and uses the Message Passing Interface (MPI) protocol for the parallel communications between the processes. Krylov subspace methods and multigrid accelerator routines are based on the Portable, Extensible Toolkit for Scientific Computation (PETSc) library[START_REF] Balay | PETSc Web page[END_REF]

Figure 3 . 1 :

 31 Figure 3.1: Numerical molecule for the Marker and Cell method. V and S denote the locations of the velocity components and strain rate tensor components.

( 3 . 1 )

 31 where the nodes for velocity components and pressure are also shown in different colors. In a structured mesh, each cell is identified univocally by a set of integers e = (i, j, k), while the six staggered nodes corresponding to the faces centers are defined by e ± e n with n = 1, 2, 3 and e 1 = (1/2, 0, 0), e 2 = (0, 1/2, 0), e 3 = (0, 0, 1/2).

3 )

 3 In eq.(3.3), φ represents the cell average value of φ. The values of the surface integrals are approximated by the midpoint rule and the values of ψ and φ are interpolated at the faces centers from the average values φ. For second order accuracy the cell-centered values of φ and its cell average value are equivalent and for this reason the overline symbol shall be dropped for the rest of this section. Eq. (3

Figure 3 . 2 : 4 r, min 3 4 r + 3 8 , 1 4 r + 3 4 3 . 1 :

 324331 Figure 3.2: Definition of the stencil for the TVD interpolation of φ.

e φ e ∆x∆y∆z ( 3 . 8 ) 9 )

 389 The evolution equation for the global fluid composition is obtained by applying the operator (3.8) to the discrete eq. (3.3): Since each cell face is counted twice with opposite normal directions, the fluxes contributions in the double summation term of eq. (3.9) cancel each other, leaving only the contribution from the boundary terms. If no fluxes are prescribed on the boundary, the numerical analog of eq. (2.7), expressing the conservation of the phase field variable, is obtained. d dt e φ e ∆x∆y∆z = o (δ) (3.10)

11 )

 11 The evolution equation for K is then obtained by scalar multiplication of the discrete momentum equation by the velocity field. The variation of kinetic energy due to the phase-field potential is thus represented by the surface tension forcing term pre-multiplied by V . dK dt surf.tens.= (V , φ∇ψ) = ∆v [

(f, g) = ∆v e f e g e ( 3 . 13 )e ψ e 3 n=1[

 3133 The rate of change of free energy is obtained by multiplying the discrete phase-field eq.(3.4) by the chemical potential. The advective term, multiplied by ψ, expresses the rate of change of free energy due to the action of the velocity field:df dt adv. = (ψ, ∇φV ) = V • nφ∆s] e±en (3.14)It is easily seen that for eq. (3.14) and eq. (3.12) to be equal and opposite it suffices that values of φ interpolated at the cell faces in the surface tension term of the momentum balance equation are consistent with numerical scheme used in the phase field equation to treat the advective fluxes. For this reason the same TVD reconstruction schemes used in eq. (3.2) must also be used for the discretization of the distributed surface tension forces.3.1.2 Treatment of the viscoelastic stressThe momentum equation also presents another delicate issue. Since the constitutive equations are solved separately in a decoupled way, the divergence of the viscoelastic part of the stress tensor is treated explicitly in the momentum equation. This may lead to a loss of ellipticity for small values of the solvent fraction parameter and cause numerical instability due to the explicit right hand side of the equation becoming dominant. For this reason, a stabilizing term is added to both sides of the momentum equation in the form of the divergence of a Newtonian stress, in a similar fashion as the Discrete Elastic/Viscous Stress Split (DEVSS) technique developed by Guénette

Figure 3 . 3 :

 33 Figure 3.3: Velocity red-black decoupling occurring for staggered arrangement of velocity and stress fields.

  (3.15) is taken as the divergence of the edge-centered strain rate tensor times the viscosity while the one on the right hand side is taken as the divergence of the cell-centered strain rate tensor times the viscosity. The parameter ̺ s represents the stabilizing viscosity and N is the non linear discrete operator arising from the discretization of the convective term of the momentum balance equation. The discrete gradient operator G p,v , defined on the space of pressure, and the discrete divergence D v,p , defined on the velocity space, are each other's transposed. This property is due to the fact that the equations are discretized on uniform cartesian meshes and it is relevant since it makes the discrete Stokes problem symmetric.

  of applying the Newton-Raphson procedure globally, each unknown is updated iteratively in a Gauss-Siedel cycle by solving a local non linear problem. Given the solution at the k th iterate the local residual and jacobian are computed as: r e = rhs e -φ k e -F e φ k j e = 1 + ∂F e (φ) ∂φ e φ k (3.22) and the local value of φ is then updated φ k+1 e = φ k e -

Figure 3 . 4 :

 34 Figure 3.4: Convergence history of the non linear GS/NGMRES method applied to the phase-field problem.

Figure 3 . 5 :

 35 Figure 3.5: Schematics of the restriction and prolongation operations on a staggered grid. Cell-centered residual in the node I is obtained by averaging fine grid values of nodes A to H (left). Staggered residual component in the node E is obtained by averaging the fine grid values corresponding to nodes A to D

  is denoted by I p n n+1 . The prolongation operator I p n+1 n is instead defined by assigning to each fine level cell the same value of the corresponding coarse level cell. The latter corresponds to a nearest-neighbor interpolation. d n i,j,k = d n-1 i/rx,j/ry,k/rz (3.29) The staggered momentum residual can be restricted by averaging each component over the corresponding cell face. For the U component, for example, the restriction operator I u n n+1 would be

  26) is bounded by the union of the circles whose centers and radii are respectively the matrix diagonal entries and the sum of the absolute values of the non-diagonal entries in each row. Since the diagonal entries of the matrix scale linearly with the viscosity of the fluids, the greater the viscosity difference between the fluids the more spaced out are the Gershgorin circles. The problem can be removed by employing the symmetric diagonal scaling technique. The original system Ax = b is replaced by the scaled system Ax = b, where A = S -1 AS -1 , x = Sx, b = S -1 b and S is the scaling matrix usually taken as a diagonal matrix whose entries are the square root of the diagonal entries of A. In this way the diagonal entries of the scaled matrix are all unitary and the Gershgorin circles become clustered in same region. The problem with the Stokes discrete operator is that its diagonal elements are zeros in the pressure block, so while the velocities unknowns are still scaled by S v = diag √ F ii , the pressures unknowns are scaled by the diagonal of the pressure Schur complement

  37) end for This procedure was first introduced by Vanka [58]. Full Vanka Smoothing (FVS) requires n cells inversions of 7 × 7 matrices. These are in general dense because of the mixed second derivatives in the velocity components appearing in reason of the non uniformity of the viscosity. For uniform viscosity the block F c is tridiagonal. For the direct solution of the local linear systems by singular values decomposition, the LAPACK [3] routine dgsev is employed. A variation of FVS consists of approximating the block F c by its diagonal. In Diagonal Vanka Smoothing or DVS, the local saddle-point problem can be solved in a very efficient way considering that, since the system is diagonally scaled, the diagonal of F c is the unit matrix. The same holds for the cell based Schur complement S c = G ′ c G c = 1. Hence the solution for the velocity and pressure corrections is obtained by: ∆p = (r p -D c r v ) /S c ∆v = r u -G c ∆p (3.38)

e i%Figure 3 . 6 :

 36 Figure 3.6: Schematic of a subdomain block. face, edge and corner ghost cells are represented only on one side for simplicity.

Figure 3 . 7 :

 37 Figure 3.7: Speed-up factor for the parallel solution of the SINKER problem described in section (4.2) for a 256 3 grid (non balanced) and a 384 3 grid (balanced).

chapitre ( 2 )

 2 , modélisant l'évolution spatiale et temporelle de l'enrobage de deux fluides, se base sur une technique de discrétisation spatiale en volumes finis appliquée à des maillages cartésiens uniformes de type MAC (maillage décalé) et caractérisée par une précision d'ordre deux. L'intégration temporelle s'effectue par une combinaison du schéma de Crank-Nicolson, utilisé pour l'équation d'évolution du champ de phase, et une formule de type BDF (backward differentiation formula) pour l'avancement en temps de V , p et τ . Seuls les termes de diffusion de phase et quantité de mouvement sont traités implicitement tandis que les autres termes sont extrapolés à partir des évaluations aux pas de temps précédents avec un schéma du type Adams-Bashforth. Ce choix permet de découpler les équations et de les résoudre séparément tout en gardant une précision du second ordre en temps. À chaque pas de temps le problème discret de Stokes est résolu par une technique itérative: la Méthode de Minimisation du Résidu Généralisée (GMRES) ou la Méthode du Résidu Conjugué Généralisée (GCR), se servant d'un préconditionneur de type multigrille pour accélérer la convergence. Il est basé sur une généralisation de la méthode de Gauss-Seidel (GS) dite de Box-relaxation, qui consiste à appliquer GS à l'ensemble des inconnus de vitesse et pression associés à une cellule et permet de contourner le problème lié à la présence d'éléments diagonaux identiquement nuls associés à la pression. L'avancement du champ de phase demande aussi méthode résolution itérative d'un système nonlinéaire par une méthode de Newton-Krylov qui offre une convergence quadratique en appliquant la méthode GMRES à un problème non-linéaire et une stratégie de préconditionnement de type Newton-Raphson locale.

Figure 4 . 1 :

 41 Figure 4.1: Comparison between analytical and numerical solution for Poiseuille flow of a Newtonian fluid in a square section. Velocity profiles are plotted at equally distanced positions along y (a). Numerical error as function of the mesh resolution showing the second order accuracy of the scheme (b).

  3.1), with a layer and flow rate repartition of

Figure 4 . 2 :

 42 Figure 4.2: Comparison between analytical and numerical solution for the Oldroyd-B flow between parallel plates with β = 0.1. (a) Velocity profile; (b) Shear stress; (c) Normal stress. (d) phase field.

Figure 4 . 3 :

 43 Figure 4.3: Error convergence for the normal stress and the phase field.

  solver. The SINKER configuration is sketched in fig. (4.4). The domain consists of a unit cubic box. If the origin is placed at the center of the box, the high viscosity central region is defined by: -0.15 ≤ x ≤ 0.15, -0.15 ≤ y ≤ 0.15, -0.15 ≤ z ≤ 0.15 (4.3)

Figure 4 . 4 :

 44 Figure 4.4: Schematic of the SINKER problem also showing streamlines for ∆η = 10 -3 . The central block in red is made of the more viscous and more dense fluid

  .1), in term of computational time and number of iterations needed for convergence. The termination criterion used to check convergence is given by ||r||/||f || < tol, where r is the residual of the system and f the right hand side. The convergence tolerance is fixed at 10 -10 . Computational times are reported with the sole purpose of comparing the different smoothers even though they are not fully representative of the real solver timing performances. Indeed, in transient simulations less stringent convergence tolerances are imposed and also the iterative solver becomes progressively more efficient as the solution is advanced in time because the initial guess solution, taken as the previous time step solution, becomes closer to the current solution. The convergence rate is estimated as [log 1 0 (||r 0 ||) -log 1 0 (||r it ||)] /it, where it is number of iterations needed for convergence.

Figure 4 . 5 :

 45 Figure 4.5: Convergence history of the SINKER problem solved with GCR+FVS on 64×64×64 grid with a viscosity jump defined over 3.5 cells and different values of the viscosity contrast.

  (4.5), the convergence history is plotted against the iteration number for the GCR+FVS combination. The difference between ∆η = 1 and ∆η = 10 -7 is only in one additional iteration required for the latter case.The effect of variable mesh size is illustrated in fig. (4.6) where the sinker problem is solved for h = 1/32, h = 1/64 and h = 1/128 and a viscosity contrast of ∆η = 10 -3 . Convergence to the fixed tolerance of 10 -6 is reached for 6, 7 and 8 iterations respectively, showing a weak sensitivity to the discretization parameter. The solution times are 3.18, 29.5 and 275.7 seconds respectively, indicating a good scalability of the time per iteration with respect to the number of unknowns.

Figure 4 . 6 :

 46 Figure 4.6: Solver sensitivity to h, for ∆η = 10 -3 and a viscosity jump defined over 3.5 cells.

Figure 4 . 7 :

 47 Figure 4.7: Schematic of the domain used for the calculation of secondary flow in square ducts.

  presented in fig. (4.7).For a direct comparison with the results of Yue, the same spatial resolution of 50×50 is used along x and y, even though convergence of the results is already observed at lower resolutions. Only five points are used to discretize the domain dimension L z because the flow is fully developed in this direction.

Figure 4 . 8 :

 48 Figure 4.8: Intensity of secondary flows in a square duct scaled by the intensity of the main flow (a), and by L/λ H (b). Comparison with numerical data of Yue et. al. [66, fig. 3]

Figure 4 . 9 : 1 SS 1 S S √ u 2 + v 2

 491122 Figure 4.9: Stress difference τ eyy -τ exx for De = 1, a = 0.5 β = 0.1. (a) Calculation by Yue et al. [66, fig. 2(c)] (b) Present simulation. (c) Streamlines of secondary flow.

Figure 4 . 10 :

 410 Figure 4.10: Schematic of the domain used for quasi-2d calculation. The aspect ratio of the rectangular section is equal to 4.

d 1 Original 10 -5 5 . 3 • 10 4 2500 1 -0 3 • 35 Rescaled 5 . 6 • 35 Table 4 . 5 :

 110531335563545 10 -3 0.1 -1 0.5 -0 0.10 -5 2.8 • 10 5 1.4 • 10 4 5.6 -0 3 • 10 -3 0.1 -1 0.5 -0 0.Flow parameters as originally reported in the work of Yue et al.[START_REF] Yue | Elastic encapsulation in bicomponent stratified flow of viscoelastic fluids[END_REF] (first row), and after rescaling (second row).

  fig.(4.11), presents central layer breakup with the formation of the lateral blobs as

Figure 4 . 11 :

 411 Figure 4.11: Temporal evolution of anomalous encapsulation (rescaled time). (a) t = 0 . (b) t = 27.85 . (c) t = 50.13 . (d) t = 83.75 . (e) t = 108.61 . (f) t = 136.46 .

Figure 4 . 12 :

 412 Figure 4.12: Comparison between the numerical solution in fig [65, fig. 10(a)] at time t = 76.21 (bottom) and the present solution at t = 136.46 (top). Temporal evolution of the interface is slower due to the secondary motions being weaker because of the greater Deborah number in our simulation.

Figure 4 . 13 :

 413 Figure 4.13: Shear viscosity material function of the Giesekus model for PS 1161 and PS 4801

04 and β 2

 2 = 0.13. The capillary number is computed on the base of an estimated surface tension of 10 -2 N/m. The Péclet number and the Cahn number are chosen on the base of purely numerical considerations: Cahn number is dictated by the fact that the interface should be at least 3.5 cells wide on the coarsest mesh, while the Péclet number is set to be 1000. Tuning the P e at a lower value would excessively damp the flow while for higher values the interface would be smeared out due to artificial diffusion.

Figure 4 . 14 :

 414 Figure 4.14: Schematic of the domain used for three-dimensional calculation of encapsulation.

L

  e and extending from z = 0 to z = L r . In this way, at the point where the two fluids physically meet, the velocity and stress profile are more naturally developed. For this second case L r = 5 and L z = 15. In both cases only half of the original domain is considered for symmetry reason and open boundary conditions are imposed for velocity and stress fields at the outlet. The flow rate ratio is unitary.

Figure 4 . 15 :

 415 Figure 4.15: Schematic of the three-layers coextrusion configuration.

Figure 4 .

 4 Figure 4.16: Mid-plane view of three layers coextrusion of Newtonian fluids with R η = 0.5 and Q 1 /Q 2 = 1 (flow rate ratio) at t = 0.25. (a) h = 0.03125; (b) h = 0.015625; (c) h = 0.0078125.

Figure 4 .

 4 Figure 4.17: Mid-plane view of three layers coextrusion of Newtonian fluids with R η = 0.1 and Q 1 /Q 2 = 5.25 (flow rate ratio) at t = 0.5. (a) h = 0.03125; (b) h = 0.015625; (c) h = 0.0078125.

Figure 4 . 18 :Figure ( 4 . 18 )

 418418 Figure 4.18: Temporal evolution of the PS 4801 layer. (a) t = 0, (b) t = 5, (c) t = 10, (d) t = 35, (e) t = 60.

( 4 .

 4 18(b)), is due to the discontinuity in the initial condition for the velocity, which is uniformly zero in the interior of the domain and finite at the inflow boundary. The instantaneous adaptation of the flow to the boundary conditions produces an abrupt expansion of the more viscous component. This disturbance propagates rapidly toward the exit of the die. The real encapsulation is a much slower mechanism taking roughly 60 non dimensional time units to reach a steady state. After this time, changes in the interface position are almost undetectable even if full encapsulation has still not been reached. This is not surprising since secondary flows are always weaker than the primary flow by one or even two orders of magnitude, as seen in section (4.3). In order to achieve full encapsulation under the same flow conditions, a die of a length equal at least to 30 should be considered. For Newtonian fluids the expansion of the more viscous component is the only observable effect since the second normal stress difference is zero and the absence of any cause for recirculation prevents the onset of encapsulation.In fig.(4.19), interface positions are compared between the two configurations considered. Since the two fluid components already have a developed stress profile, encapsulation develops more rapidly in case the separation plate is included. Apart from the bulge appearing at the centerline of the section, which might also be due to the unphysical inflow conditions, the profiles are similar in the two cases except that they are shifted along the axis z. Experimental data (courtesy of Teixeira, 1996[START_REF] Teixeira | Étude expérimentale et numérique du phénomène d'enrobage en coextrusion[END_REF]) are reported in fig. (4.20(a)). The first row of slices corresponds to the interface profiles at positions from z = 0 to z = 7.38. The encapsulation profile obtained by numerical simulation is compared with the experimental data in fig.

Figure 4 . 19 :

 419 Figure 4.19: Interface profiles for configurations with (right) and without (left) the separation plate. (a,e) z = 0.5, (b,f) z = 2.5, (c,g) z = 5, (d,h) z = 8

Figure 4 . 54 Résumé du chapitre 4

 4544 Figure 4.20: (a) Experimental interface profiles from z = 0 to z = 24. Data courtesy of Teixeira, 1996 [44] (b) Comparison between the numerical encapsulation at z = 5 and the experimental data at z = 5.54

Figure 5 . 1 :

 51 Figure 5.1: Schematic side view of the computational domain used for the simulations of the first regime.

Figure 5 . 2 :Figure 5 . 3 :

 5253 Figure 5.2: Comparison of two-layers coextrusion flows: η 2 /η 1 = 0.2 (left), η 2 /η 1 = 0.05 (right). t = 0 (a,e). t = 2 (b,f). t = 4 (c,g). t = 6 (d,h).

Figure 5 . 4 :

 54 Figure 5.4: Interface profile in the cross section plane z = 1 at t = 6. η 2 /η 1 = 0.2 (a). η 2 /η 1 = 0.05 (b).

Figure 5 . 5 :

 55 Figure 5.5: Interface profiles in the symmetry plane for different values of the viscosity ratio.

Figure 5 . 6 :

 56 Figure 5.6: Interface position and entry length as functions of the viscosity ratio (a). Dimensionless pressure gradient and viscous dissipation as functions of the viscosity ratio (b).

Figure 5 . 7 :

 57 Figure 5.7: Evolution of the interface shape due to a sinusoidal flow rate perturbation in a two-layers stratified flow with Re = 10 -6 , η 2 /η 1 = 0.1. t = 4 (a). t = 5.5 (b). t = 7 (c). t = 8.5 (d). t = 10 (e). t = 11.5 (f)

Figure 5 . 8 :

 58 Figure 5.8: Lateral view of the computational domain used for the simulations of the second regime of encapsulation (a). Local coordinate system for the definition of N 2 at the wall (b)

2 . 2 )

 22 . Oldroyd-B model fails in the prediction of a second regime of encapsulation since it exhibits no second normal stress difference and the same consideration made for the Newtonian model can be drawn for Oldroyd-B fluids. To be able to observe the second regime, the Giesekus model is adopted in this section. With this model the second normal stress difference is controlled by the mobility coefficient a which corresponds to the ratio -N 2 /N 1 in the limit of zero shear rate. A known issue of the Oldroyd-B model, and by extension of the Giesekus model as a tends to zero, is the presence of an unbounded extensional viscosity. When imposing a boundary inflow condition like the one used in section (5.2.1), a double Poiseuille profile, the fluids are subject to a severe elongation in the region where the two streams merge, since the point of the boundary where the two velocity profiles meet is characterized by a virtually infinite ∂w ∂z . In numerical simulations the value ∂w ∂z is finite but

  (5.8(b)), at the wall, where the flow direction is z ≡ 1, the velocity gradient direction is x ≡ 2 and the third neutral direction is y ≡ 3, one has N 2 = τ xx -τ yy . Similarly, at the symmetry plane z ∼ = 1 , y ∼ = 2 and x ∼ = 3, and then N 2 ≈ -(τ xx -τ yy ).

Case 2 : η 2 /η 01 = 2 Figure 5 . 9 :Figure 5 . 10 :Figure 5 . 11 :

 22259510511 Figure 5.9: Shear viscosity as a function of the shear rate for Newtonian and Giesekus fluids.

134 Chapter 5 .Figure 5 . 12 :

 5512 Figure 5.12: Encapsulation level for the case N EW T GIE1 at W i = 0.2154 (a), as defined in the schematic (b).

of η 0 G

 0 /η N = 0.5 for which the viscosity of the Newtonian fluid is always greater than the one of the Giesekus fluid. As for the previous case, the secondary flow intensity, shown in fig. (5.13(a)), attains its maximum around W i = 0.21 and the analysis of the secondary flow fields in fig. (5.13(b)), (5.13(c)) and (5.13(d)), indicates that the Newtonian fluid encapsulates the Giesekus one in any case. This is striking in the three-dimensional visualization of the Newtonian layer in fig. (5.14(c)) to (5.14(k)).

5 and a 1 , 2 = 0. 6 .

 126 This choice produces a crossover of the second normal stress difference coefficients of the two fluids and a double crossover in the viscosity, as shown in fig. (5.15(a)) and (5.15(b)).

Figure 5 . 13 :

 513 Figure 5.13: Intensity of the secondary flow as a function of the Weissemberg number for the configuration N EW T GIE2 (a). Visualization of the secondary flow in the left half of the square cross section. Newtonian fluid is red (upper) and Giesekus fluid is blue (lower). W i = 0.0278 (b). W i = 0.2154 (c). W i = 1.6681 (d).

Figure 5 . 14 :Figure 5 . 15 :

 514515 Figure 5.14: Case N EW T GIE2 at W i = 0.2154. Inlet velocity profile at indicated x sections (a) and normal stress (b). Visualization of the deformation of the Newtonian fluid layer. t = 0 (c). t = 5.25 (d). t = 10.5 (e). t = 15.75 (f). t = 21 (g). t = 26.25 (h). t = 31.5 (i). t = 37.75 (k)

1 . 2 .Figure 5 . 16 :

 12516 Figure 5.16: Intensity of the secondary flow as a function of the Weissemberg number for the configuration GIEGIE2 (a). Visualization of the secondary flow in the left half of the square cross section. Fluid 1 is in red (upper) and fluid 2 fluid is in blue (lower). W i = 0.2154 (b). W i = 1.1 (c). W i = 3 (d).

A negative N 2 Figure 5 . 17 :

 2517 Figure 5.17: Case GIEGIE2 at W i = 3. Inlet velocity profile at indicated x sections (a) and normal stress (b). Visualization of the deformation of the fluid 2 layer. t = 0 (c). t = 7.5 (d). t = 15 (e). t = 22.5 (f). t = 30 (g). t = 37.5 (h). t = 45 (i). t = 52.5 (j). t = 60 (k).

Figure 5 . 18 :

 518 Figure 5.18: Case GIEGIE2 at W i = 1.1. Inlet velocity profile at indicated x sections (a) and normal stress (b). Visualization of the deformation of the fluid 2 layer. t = 0 (c). t = 7.5 (d). t = 15 (e). t = 22.5 (f). t = 30 (g). t = 37.5 (h). t = 45 (i). t = 52.5 (j). t = 60 (k).

Figure 5 . 19 :

 519 Figure 5.19: Schematic of mechanism driving the interface deformation at the wall for the cases N EW T GIE1 and N EW T GIE2.

Figure 5 . 20 :

 520 Figure 5.20: Schematic of the three-layers configuration.

h 0 [Figure 5 . 21 :

 0521 Figure 5.21: Three-layers configuration: comparison of the interface position in the symmetry plane for different viscosity ratios.

  2.3) indicating that the Newtonian fluid always encapsulates the Giesekus fluid. In this way N 2 is always greater in one fluid irrespective of the W i of the flow, and the effect of geometry can be clearly evinced since the way the jump of N 2 influences the interface deformation for this fluids combination has already been elucidated. All the flow conditions and simulation parameters are the same as the one of the corresponding two-layers flow, except for the presence of a double symmetry boundary condition allowing the restriction of the domain to only one quarter of the original section. As for the previous cases, the quasi2D flow with the fixed interface (corresponding to the inflow boundary condition) is inspected first in order to verify how predictive is the secondary flow pattern with respect to the actual three-dimensional deformation of the interface. Moreover, this step allows the determination of the W i yielding the strongest recirculation and thus the fastest encapsulation rate.From the analysis of the secondary flow intensity curve plotted against the Weissemberg number in fig.(5.23(a)) it appears that when the core layer is made of the Giesekus fluid (configuration B -A -B) the curve is shifted forwards and the maximum of the intensity occurs at a higher W i than for the opposite case when the core layer is made of the Newtonian fluid (configuration A -B -A). This is far more evident for the rectangular section. This effect is of difficult interpretation

Figure 5 . 22 :Figure 5 . 23 :

 522523 Figure 5.22: Comparison of three-layers coextrusion flows (only half channel): η 2 /η 1 = 0.1 (left), η 2 /η 1 = 10 (right). t = 0.75 (a,e). t = 1.5 (b,f). t =

( 5 .

 5 23(b)) is that the encapsulation characteristic time, estimated from the average secondary flow intensity, is far greater than the polymer natural relaxation time, meaning that the elasticity effects are negligible for this case too.

Fig. ( 5 Figure 5 . 24 :

 5524 Fig. (5.24) shows the secondary flow pattern in the square section and confirms the observations made for the two-layers flow. For both configurations the flow is directed from the Newtonian fluid to the Giesekus fluid near the wall and from the Giesekus fluid to the Newtonian fluid near the symmetry plane. The threedimensional results confirm this tendency as shown by fig. (5.26) and (5.27).

Figure 5 . 25 :

 525 Figure 5.25: Secondary flow for the three-layers configuration in the rectangular section at the Weissember numbers yielding the maximum recirculation intensity. A -B -A and W i = 0.2154 (a). B -A -B and W i = 0.5995 (b). B is the Newtonian fluid and A is the Giesekus fluid.

Fig. ( 5

 5 Fig. (5.29) and (5.27) show how neither in the rectangular nor in the square section B -A -B configurations the contact lines come in contact within the first 60

Figure 5 . 26 :

 526 Figure 5.26: Interface deformation for the three-layers configuration A-B-A with square section at W i = 0.2154. t = 7.5 (a). t = 15 (b). t = 22.5 (c). t = 30 (d). t

Figure 5 . 27 :

 527 Figure 5.27: Interface deformation for the three-layers configuration B-A-B with square section at W i = 0.2154. t = 7.5 (a). t = 15 (b). t = 22.5 (c). t = 30 (d). t = 37.5 (e). t = 45 (f). t = 52.5 (g). t = 60 (h).

Figure 5 . 28 :

 528 Figure 5.28: Interface deformation for the three-layers configuration A -B -A with rectangular section, L x /L y = 2, at W i = 0.5995. t = 7.5 (a). t = 15 (b). t

  t

Figure 5 . 29 :

 529 Figure 5.29: Interface deformation for the three-layers configuration B -A -B with rectangular section, L x /L y = 2, at W i = 0.5995. t = 7.5 (a). t = 15 (b). t

  La méthode du champ de phase est une alternative très intéressante aux méthodes classiques de suivi d'interface pour la simulation numérique des écoulements en coextrusion en raison de sa prédisposition à gérer à la fois le mouvement du point triple de contact d'une interface avec une paroi solide et le changement de topologie des sous-domaines occupés par chaque fluide causé par l'enrobage. Malgré que le code numérique implémenté ne représente qu'une version préliminaire et nécessite des développements ultérieurs afin de pouvoir atteindre les performance requises pour l'obtention de résultats quantitativement significatifs, la méthode proposée s'est révélée capable de modéliser les principaux facteurs influençant le phénomène de l'enrobage. Par l'analyse des solutions 3D obtenues pour les cas de coextrusion de deux ou trois couches dans un filière à section carrée ou rectangulaire, nous sommes parvenus à caractériser la nature des deux régimes d'enrobage: le premier régime comporte une réadaptation des épaisseurs des couches provoquée par le déséquilibre des gradients de pression sur l'interface, ce qui est directement lié aux rapports de viscosité entre les fluides et aux débits; alors que le deuxième régime est dû aux différences de seconde différence de contraintes normales entre deux couches adjacentes. Le fluide présentant le N 2 plus important en valeur absolue présente localement une interface convexe. Une vraie limite rencontrée dans le modèle du champ de phase se trouve dans le fait que le vrai "moteur" de l'enrobage -le saut de N 2 localisé sur l'interface-souffre de la précision relativement faible du modèle dans cette région de l'écoulement. Deux possibilités sont envisagées pour régler ce problème. Une première solution serait de modifier l'équation de Cahn-Hilliard en rajoutant des termes spécifiquement conçus pour contraster ces effets non-désirés, l'autre est d'utiliser un Raffinement Adaptif de Maillage (AMR) permettant d'augmenter la résolution locale du maillage et donc la précision de la solution uniquement près de l'interface. Cette amélioration pourrait permettre l'extension de l'étude de l'enrobage aux points suivants: 1. influence des différentes conditions aux limites sur le comportement du point triple à paroi, 2. influence des effets de tension de surface, 3. effet d'une géométrie variable, 4. effet d'un modèle viscoélastique de type multi-mode sur le degré d'enrobage. Ce dernier point serait particulièrement intéressant car la présence d'un maximum absolu dans la courbe d'intensité des écoulements secondaires en fonction du nombre de Weissemberg W i est associé au temps de relaxation du polymère. Donc, pour un modèle à un seul temps de relaxation il n'existe qu'une seule condition qui maximise l'intensité de l'écoulement secondaire et donc le degré d'enrobage, alors que pour un modèle multi-modes il pourrait y avoir un intervalle de W i dans lequel le phénomène d'enrobage est plus présent.

  e 1 ; e i + e 1 ; e i -e 2 ; e i + e 2 ; e i -e 3 ; e i + e 3 ]

% Define the subset of the pressure unknowns surrounding the cell e i s p i = [e i -2e 1 ; e i + 2e 1 ; e i -2e 2 ; e i + 2e 2 ; e i -2e 3 ; e i + 2e 3 ] % Restrict the discrete diffusion gradient and divergence operators

Table 4 .

 4

	∆t Re h 2 ≈ 10 31						Smoother
	∆η		DVS			FVS			DGS			None
		it.	time	C.R	it.	time	C.R.	it.	time	C.R.	it.	time	C.R.
	10 0	11 25.8 0.843 11 43.3 0.896 12 26.6 0.778 33 33.2 0.297
	10 -1	11 29.7 0.839 11 48 0.889 13 29.7 0.73 38 40.5 0.257
	10 -2	12 32.4 0.812 11 48.1 0.857 22 50.5 0.437 44 46.9 0.225
	10 -3	12 32.3 0.772 12 52.4 0.834 34 79.7 0.289 -	-	-
	10 -4	13 35.3 0.729 12 52.5 0.823 -	-	-	-	-	-
	10 -5	14 37.8 0.694 12 52.5 0.818 -	-	-	-	-	-
	10 -6	15 42.7 0.655 12 52.5 0.81	-	-	-	-	-	-
	10 -7	15 42.7 0.649 12 52.5 0.805 -	-	-	-	-	-

1: Convergence properties of the three proposed smoothers: Iteration count (it.); elapsed time (time) and convergence rate (C.R). Solutions are computed for the sinker problem with h = 1/64 and w t /h = 3.5 and varying the viscosity ratio.

Table 4 . 2

 42

	∆t Re h 2 ≈ 10 31					Smoother
	w t /h		DVS			FVS			DGS
		it.	time	C.R	it.	time	C.R.	it.	time	C.R.
	2	28 70.2 0.347 33 133 0.298 -	-	-
	3.5	13 35.3 0.729 12 52.5 0.823 -	-	-
	5	12 32.7 0.8 11 48.6 0.857 24 55.3 0.415
	10	11 30 0.839 11 48.9 0.891 13 31.5 0.761

very efficient for high Re values as reported in table

(4.3)

. : Effect of interface thickness. Iteration count (it.); elapsed time (time) and convergence rate (C.R). Solutions are computed for the sinker problem with h = 1/64 and ∆η = 10 -4 .

Table 4 .

 4 4: Effect of interface thickness. Iteration count (it.); elapsed time (time) and

convergence rate (C.R). Solutions are computed for the sinker problem with ∆x = ∆y = h = 1/32, ∆η = 10 -1 , w t /h = 3.5,and varying the streching factor ∆z/h.

Table 4 .

 4 6: Giesekus model parameters for PS 1161 and PS 4801 Material properties and flow parameters The two polymers considered are two polystyrenes (PS) under the reference names of Lacqrene 4801 and Lacqrene 1161. PS 4801 is less viscous than PS 1161 as depicted at fig. (4.13). The experimental data are fitted by the shear viscosity material function of the Giesekus model.

	Polymer η p [P a s] η s [P a s]	λ [s]	a
	PS 1161 6128.66	254.53	0.8401 0.5246
	PS 4801 1235.85	182.202 0.46571 0.6910

Résumé du chapitre 1 L'objectif principal de cette thèse est de comprendre le mécanisme générant le phénomène de l'enrobage et, en particulier, de séparer les facteurs influençant le premier régime, consistant en une adaptation des différentes épaisseurs des couches aux conditions de gradient de pression en entrée de filière, de ceux qui déterminent le deuxième régime, c'est-à-dire le véritable enrobage. La littérature suggère une interprétation de ce phénomène essentiellement liée à la présence d'une différence de viscosité des fluides qui favoriserait la tendance du fluide le moins visqueux à migrer vers les zones de l'écoulement présentant un plus fort taux de cisaillement, de manière à minimiser la dissipation d'énergie cinétique due à la viscosité. Toutefois, la physique de la transition d'une configuration stratifiée à une enrobée reste toujours à clarifier, car un simple critère basé sur le rapport de viscosité ne suffit pas à expliquer certaines observations expérimentales comme par exemple la totale absence d'enrobage pour les fluides Newtoniens ou bien l'existence de conditions d'écoulement pour lesquels le fluide plus visqueux enrobe le moins visqueux. La méthodologie numérique adoptée pour l'analyse des écoulements de coextrusion est inspirée des travaux de Mompean at. al [41] pour la conception d'un code de calcul basé sur la méthode des volumes finis appliquée aux maillages cartésiens uniformes. La présence de deux fluides est prise en compte par le modèle du champ de phase, remplaçant l'interface mince existant entre deux fluides par une interface diffuse. Ce modèle a été utilisé avec succès par Yue et al. en [65] pour simuler l'enrobage dans l'hypothèse d'écoulement quasi2d avec la méthode des éléments finis couplée avec une technique de raffinement de maillage adaptatif. L'élément innovant par rapport à cette dernière solution est représenté par l'extension de la méthode au cas à trois dimensions. Ce passage demande une attention particulière sur certains aspects de la modélisation comme par exemple l'implémentation d'un solveur de

Acknowledgements

The author of this thesis would like to acknowledge the funding program CIBLE of the region Rhône-Alpes for financing this work; the Centre Informatique National

List of Tables

Validation of the code

The following chapter is dedicated to the validation of the numerical code for the simulation of viscoelastic two-phase flows.

Accuracy testing

The second order space accuracy of the finite volume scheme adopted is verified against analytical solutions available for simple flow configurations. The first case is represented by a stationary fully developed flow in a rectangular duct whose section is defined as [-L x /2, L x /2] × [0, L y ]. For Newtonian fluids the solution for the axial velocity w, as reported in [START_REF] Mortensen | Reexamination of hagenpoiseuille flow: Shape dependence of the hydraulic resistance in microchannels[END_REF], is:

while u and v are nil. In the expression (4.1), ∆p L and η denote respectively the pressure gradient and the dynamic viscosity. The simulations are performed on a very long duct (10 3 L x ) in order to allow the full development of the flow, imposing a constant inflow velocity profile and inspecting the solution at the outflow section.

The resolution in the cross section plane is taken as h = 1/2 n with n = 5, 6, 7, 8, while it is kept constant at 1/16 along z. The comparison between the analytical and numerical solutions is shown in figures (4.1(a)). The error is computed as the L 2 norm of the difference between the numerical and analytical solutions. The convergence rate is estimated as log 2 (err n ) -log 2 (err n+1 ). The asymptotic convergence rate calculated for the two cases is approximately 1.99, as shown in fig. region as observed in all the numerical simulations. The same considerations can be extended for the interface point in the midplane, where τ xx -τ yy corresponds to -N 2 and the interface normal lies in the third direction. In this case the stronger the N 2 the greater the compressive stress, hence the net force on the interface pulls it in the region of lower N 2 . In contrast with what is observed for the contact line movement, the flow direction in the midplane is not always consistent with the sign of the N 2 . This is a consequence of the quasi2D assumption. Since all the derivatives in the z direction are zero, the streamlines in the cross section plane are constrained to form closed loops by the incompressibility condition. This implies that the net flow rate across the initial flat interface must be zero:

The sign of v results from the competition of two factors: the local sign of the normal stress jump and the global constraint imposed by the incompressibility. At the midplane, the N 2 jump is usually weaker than at the wall so the flow is most influenced by this latter condition. This analysis explains the discrepancy observed of the flow simulation GIEGIE2 at W i = 3 between the local flow direction observed at the inflow section and the actual three-dimensional deformation of the interface in the opposite direction. Indeed, in three-dimensional flow the fluid stream is allowed to accelerate in the z direction, so the condition expressed by the relation (5.4) does not hold anymore. This consideration poses another limit to the validity of quasi2D

simulations which are able to correctly predict the interface movement at the wall but may fail to describe the interface deformation in the interior domain.

Three-layers configuration

When more than two layers are coextruded, the effect of geometry comes into play and the way the different layers are alternated has a role in the determination of interface final state too. A simple way to illustrate this effect is offered by the two-dimensional three-layers channel flow. As depicted in fig. (5.20), showing a half of the channel, the two layers have an initial repartition and flow rate ratio of