N Lebaal
email: lebaal@enstimac.fr

F M Schmidt

S Puissant
email: stephan.puissant@insic.fr

O ptimization of 3D die extrusion using response surface method

Keywords: Polymer extrusion, Response surface method, Optimization, Kriging Interpolation, REM3D

The primary objective of the geometrical design of extrusion dies in polymer processing is to obtain a uniform velocity distribution across the die exit. A design procedure for complex coat-hanger die is presented. While optimizing the exit velocity distribution, geometric constraints are applied. Three dimensional extrusion simulation software REM3D ® is used to simulate the flow in this flat die. An objective function is defined as the global relative between velocity in exit die and the average exit velocity. This objective function is minimized by varying the flow channel cross-section. For this minimization the global response surface method with Kriging interpolation is used.

INTRODUCTION

The design of dies for polymer extrusion often involves trial and error corrections to the die to achieve uniform flow at the exit. Manual correction to die geometry is a time consuming and a costly procedure. A flat die (Fig. 1) is commonly used to extrude thermoplastics thin sheets [START_REF] Sun | Effect of elongational viscosity on the flow in flat die[END_REF][START_REF] Michaeli | Approch of automatic extrusion die optimisation[END_REF][START_REF] Reifschneider | Automated sheet Die design[END_REF]. If the channel geometry in a flat die is not designed properly, the velocity at the exit of the flat die may not be uniform [START_REF] Chen | Optimization of the Coathanger manifold via computer simulation and orthogonal array method[END_REF][START_REF] Wang | The flow distribution of molten polymers in slit dies and coathanger die through three-dimensional finite element analysis[END_REF]. A non-uniform velocity distribution at the die exit may lead to a variation in the sheet thickness across the width of the die. Since a tight control on thickness is required for a high quality plastic sheet, the ultimate goal of this work is to optimize the die channel geometry in a way that a uniform velocity distribution is obtained at the die exit [START_REF] Reddy | Design and optimisation of three dimensional extrusion dies using adaptative finite element method[END_REF][START_REF] Nóbrega | Flow Balancing in extrusion dies for thermoplastic profiles[END_REF][START_REF] Nóbrega | Sensitivity of flow distribution and flow patterns in profile extrusion dies[END_REF]. Prior works in sheet die optimization have involved the use of lubrification approximations of the momentum equations [START_REF] Yu | A simple numerical approach for the optimal design of an extrusion die[END_REF][START_REF] Smith | An optimisation-based approach to compute sheeting die designs for multiple operating conditions[END_REF][START_REF] Xiaorong | Optimal design for polymer sheeting dies[END_REF][START_REF] Ettinger | Parameterization and optimisation strategies for the automated design of U PVC profile extrusion dies[END_REF]. If the geometry is more complex, a flow channel can be approximated with simple geometric sections [START_REF] Puissant | Two dimensional multilayer coextrusion flow in a flat Coat-Hanger die[END_REF]. Others have used three-dimentional analysis to design die extrusion [START_REF] Michaeli | Computer aided optimisation of extrusion dies[END_REF][START_REF] Nóbrega | Flow balancing of profile extrusion dies[END_REF] but they did not take into account the thermal dependence. Network algorithms have been developed to optimize die designs [START_REF] Hurez | A New Design procedure for profile die[END_REF] but they are difficult to apply to arbitrary shapes. Michaeli et al [START_REF] Michaeli | Approch of automatic extrusion die optimisation[END_REF] have used a combination of finite-elementanalysis and flow analysis network. To optimize the die geometry they used respectively the evolution strategy algorithm and network theory with isothermal flow. But this geometry is varied manually and the optimization by evolution strategy makes a lot of time to converge. Sun et al [START_REF] Sun | Optimization of a flat die geometry[END_REF], Smith [START_REF] Smith | An optimisation-based approach to compute sheeting die designs for multiple operating conditions[END_REF] have used respectively BFGS and SQP algorithms to optimize die geometry but the used algorithms are time consuming and it is not easy to apply them with different software and to change different geometries. Reddy et al [START_REF] Reddy | Design and optimisation of three dimensional extrusion dies using adaptative finite element method[END_REF] have used response surface method with polynomial approximation, but this approximation is not accurate and needs a lot of iteration to converge. The complexity of the polymer rheology further increases the difficulty of the die optimization problem [START_REF] Yu | A simple numerical approach for the optimal design of an extrusion die[END_REF]. If the polymer rheological behaviour is not accounted accurately while optimizing the die, the computed velocity, pressure and temperature fields are expected to have large errors. There are many works in the literature that used rheological behaviours of power law [START_REF] Yu | A simple numerical approach for the optimal design of an extrusion die[END_REF][START_REF] Smith | An optimisation-based approach to compute sheeting die designs for multiple operating conditions[END_REF][START_REF] Ettinger | Parameterization and optimisation strategies for the automated design of U PVC profile extrusion dies[END_REF][START_REF] Michaeli | Computer aided optimisation of extrusion dies[END_REF][START_REF] Hurez | A New Design procedure for profile die[END_REF]; however, this does not permit to present accurately the rheological behaviour. Sun et al [START_REF] Sun | Optimization of a flat die geometry[END_REF] have optimized a flat die and they have taken into account the elongation effect. They obtained optimal distribution but the pressure drop could not be decreased. In the present work we have developed an automatic optimization algorithm based on response surface method together with Kriging interpolation. We used REM3D ® software to compute 3D simulation of the flow in extrusion dies. This software takes into account strain rate and temperature dependence. The optimal design procedure is applied to a coat hanger die used to extrude thermoplastic sheets. The viscosity models and the values of the various parameters of these models are summarized in the next section for the polymer used in the present work. In this paper, a geometry of a flat die is optimized for an Acrylonitrile Butadiene Styrene (ABS, Astalac EPC 10000). The polymer flow is supposed purely viscous (the viscoelastic behaviour is not taken into account). According to Sun and Gupta [START_REF] Sun | An analysis of the effect of elongational viscosity on the flow in a flat die[END_REF][START_REF] Sun | Optimization of a flat die geometry[END_REF] the elongation effect does not influence the exit velocity distribution at the die exit. The simulation of extrusion is carried out on a 3D computation software by finite elements REM3D ® [START_REF] Pichelin | Finite element solution of the 3D mold filling problem for viscous incompressible fluid[END_REF]. The behaviour laws used in Rem3D give an expression of the viscosity in function of the shear rate and temperature. The thermal and rheological parameters of the ABS are given in Table 1. Carreau Yasuda/WLF viscosity model is used to characterize the temperature and shear rate dependence of viscosity [START_REF] Balasubrahman | Thermal control of melt flow in cylindrical geometries[END_REF]. It is written as:

( ) ( ) ( ) ( ) α α τ γ η η η 1 0 2 1 2 1 0 1 exp . - ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ - + - - - + - = m s s s s ref s ref T T A T T A T T A T T A & (1)
The thermal conductivity (K), density (ρ) and heat capacity (C p ) of the polymer used for the flow simulation were assumed to be constant throughout the range of temperature in the flat die. In order to save time at the optimization process and to be able to control effectively and easily the preprocessing REM3D (GLPre), the solver REM3D and design routine, it is necessary to use all the possibilities of automation offered by this code (via Matlab). Matlab makes it possible to launch the executable files. This fact enables us to couple Matlab with REM3D. Fig. 2, shows a diagram of the principle of operation and possible interactions of Matlab with REM3D/GLpre, REM3D/Solver and design routine. It is also possible to automate all the tasks usually carried out through graphic interface from the creation of FEM until the recovery of the results.

DESIGN VARIABLES AND OBJECTIVE FUNCTIONS

The variable of optimization is the depth of the channel of distribution H. Initially, it is fixed at 36.5 mm, but during the processes of optimization it is limited between 30 and 110 mm (fig. 1). It was mentioned earlier that the goal in design of a flat die is to minimize the velocity variation across the die exit. The objective function is:

∑ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = N i i v v v N J 1 2 1 ( 2
)
Where N is the total number of nodes at the die exit in the middle plane, i is the velocity at an exit node, v v is the average exit velocity.

OPTIMISATION PROCEDURE

The method of response surface [START_REF] Lebaal | Rheological parameters identification using in-situ experimental data of a flat die extrusion[END_REF] consists in the

Definition of die channel geometry

Solver REM3D Simulation of no-isothermal 3D Optimization routine Evaluation of the geometry construction of an approximate expression of objective function starting from a limited number of evaluations of the real function. The main idea is to approximate the objective function through a response surface. In order to obtain a good approximation, we used a Kriging interpolation [START_REF] Trochu | Nonlinear modelling of hysteretic material laws by dual kriging and Application[END_REF]. In this method, the approximation is computed by using the evaluation points by composite design of experiments, the number of evaluation is 5 points which permits to give a good interpolation. After the interpolation of the objective function we minimize it using SQP sequential quadratic programming method. In addition, to avoid to obtain a local minimum, we change the initial points of the SQP method by all the points of the design of experiment because the evaluation of the interpolation function does not cause problems and does not take much time. After we obtain the best minimum of this interpolation we use the weight function of Gaussian type which allows to slightly change the interpolation and gives more importance to the points which are closer to the minimum and less importance to the other points. The iterative procedure stops when the successive points are superposed with a tolerance ε=10 -3 .

RESULTS AND DISCUSSION

ow conditions, the optimal value for H is 94.9 mm. After 2 iterations, the objective function is varying quite fast at the beginning of the optimization until the 2 nd step. Its final value is about 5.3 10 -2 . However the global relative error between exit velocities in initial and optimal die, Fig. 3, is minimized for 94.7%. This optimization run represents 11 hours 23 min of CPU time on a computer Pentium IV, 3 GHz, 1Go RAM.. Given the fixed geometry constraints of the die, the rheological parameters of the ABS and the fl Velocity distributions at the exit of the optimized die in Fig. 4-B are more uniform than the exit velocity in the initial die design in Fig. 4-A, it is evident that in the initial die design the exit velocity at the border of the die is significantly higher than the velocity in the centre of the die. In Fig. 4-B, the exit velocity distributions for the optimized die design are much more uniform across the complete width of the half die. Only in the small regions near the border of this die, the exit velocity is reduced to zero to satisfy the no slip condition at the walls. The pressure drop in the initial die design is 4.423 MPa Fig. 5-A. It should be noted that the pressure drop decreases during the optimization iteration. The total pressure drop in the optimal die was decreased by 49.9% or 2.215 MPa, Fig. 5-B. With the optimized geometry, we can obtain the same pressure as with the initial geometry by increasing the flow rate by 140%. To summarize, it is possible to increase the productivity of more than 140% with the same pressure as for the initial die. To improve the uniformity of die exit velocity, the width of the channel cross-section is typically increased from 36 mm to 94.9 mm. Bad temperature distribution is observed in the initial die Fig. 6-A. However, it is noted that more polymer is flowing towards the edge because of the shear heating. Hence, the temperature is lower near the middle of the die and higher near the edge. By comparison, the shear heating in the optimal die and the temperature distribution are more uniform Fig. 6-B, because of the homogeneous distribution of the exit velocity.

CONCLUSION AND PROSPECTS

We conclude that we obtain good results using the response surface method with a very fast convergence (2 iterations), this is due to the good interpolation obtained with Kriging. The weight function and the various evaluations of the interpolation function (which does not take much of CPU time) permit to give a precise global minimum. A program for optimizing the geometry of a flat die for polymer sheet extrusion has been developed. To save time and to control easily the optimisation procedure, we have automised this program with REM3D software and design routine which we have also developed. The optimization program used the REM3D software to simulate the polymeric flow in the flat die. This program optimized successfully the 3D geometry of a flat die such that a uniform velocity and temperature distribution were obtained at the die exit with decreasing the pressure drop in the die. In our perspectives, we envisage to optimise the temperature of regulation instead of modifying the geometry while keeping the same objective function, that is to say homogenizing the exit velocity in the die.

Fig. 1 .

 1 Fig. 1. Geometry of the Flat die

Fig. 2 .

 2 Fig. 2. Optimization algorithm.

Fig. 4 .

 4 Fig.3. Error between th optimal and initial exit velocity.e average,

Fig. 5 .

 5 Fig. 5. Pressure distributions in initial and optimized die

Fig. 6 .

 6 Fig. 6. Temperature distributions in initial and optimized die.


DETERMINATION OF THE ELONGATIONAL PROPERTIESOF POLYMERS USING A MIXED NUMERICAL-EXPERIMENTAL METHOD

L. Robert, S. Hmida-Maamar, V. Velay and F. Schmidt Research Centre on Tools, Materials and Processes (CROMeP) Ecole des Mines d'Albi-Carmaux, 81013 ALBI Cedex 9, France {laurent.robert,hmidahmi,vincent.velay,fabrice.schmidt}@enstimac.fr

The biaxial properties of materials such as rubbers or polymers are often difficult to identify, and generally a simple isotropic behaviour derived from classical 1D tensile test is considered. However, biaxial properties are useful for the simulation of plastic-processing operations such as blow moulding [START_REF] Rosato | Blow molding handbook[END_REF] or thermoforming. As previously reported, the rheological behaviour and the mechanical properties of rubbers and polymers can be obtained by using a bubble inflation rheometer [2-4], or multi-axial tensile test [5]. In this work, experimental data provided by optical measurements on tensile tests and bubble inflation tests are coupled with Finite Element Method simulations for identifying the rheological behaviour. This work is actually based on natural rubber and will be extend to thermoplastic (PP and PET) materials.

A bubble inflation rheometer has been developed in the laboratory. It allows to blow under r r controlled pressure rubber or membranes [2]: a circular membrane, clamped at the rim, is inflated by applying air pressure to its bottom face (see Fig. 1, left). In the case of thermoplastic, a heating step is necessary before applying the pressure. The heating can be performed by air convection, by conduction (heating cartridge at the rim) and by IR radiation. f To perform DISC directly on the rubber bubble, several difficulties are to be solved (large level of deformation, semi-transparent aspect of the materials, t lighting, etc.).

In addition, tensile tests have been performed using DISC due to the high level of strain. Tensile test performed on standard specimen give, on the one hand, the stress/strain curve. On the The Mooney-Rivling hyperelastic behaviour has been chosen for its simplicity (2 parameters) and because it is representative of natural rubbers behaviour, at least for reasonab t le stretch ratio (no strain hardening). Fitting of the 1D stress/strain curve gives a set of parameters using eq. ( 1), numerical values obtained are C 10 = 0.136 MPa and C 01 C C = 0.097 MPa.

(

) 1
FIGURE 2. Experimental and numerical shape contours after optimization.

Then identification of the Mooney Rivling coefficients based on experimental shape contours using a Finite Element Model Updated (FEMU) procedure is performed using Abaqus ® . The cost function to minimise is defined as the norm (least squares sense) between experimental and numerical bubble shapes, and a Sequential Quadratic Programming (SQP) optimisation algorithm has been used. Results show that the Mooney-Rivling behaviour seems to be correct for small stretch ratio (< 200%) (see Fig. 2).

The second identification approach is based on a FEMU procedure using full field measure-d d ments by DISC (Vic-3D ® ) during the open hole test. Boundary conditions, cost function and preliminary results of the identification will be presented. Parameters values resulting from the identification based on the heterogeneous displacement field will be compared with those identified from the contour extractions and discussed.

FIGURE 1 .

 1 FIGURE 1. Bubble inflation rheometer (detail, left) and shape contour extraction (right).

L.

  Robert et al. other hand, open hole tensile tests are performed on drilled specimen, putting in evidence non-d homogeneous strain fields measured by the optical full field method. d

  2. Reuge, N., Schmidt, F.M., Lemaoult, Y., Rachik, M. and Abbe, F., Polym. Eng. Sci., Vol. Shaw, J. and Wineman, A., Exp. Mech. vol 46-5, 579-587, 2006. 5. Chevalier, L. and Marco, Y., Polym. Eng. Sci., vol. 42-2, 280-298, 2002. 6. Garcia, D., Orteu, J.-J. and Penazzi, L., J. Mat. Proc. Tech., vol. 125, 736-742, 2002. Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org)

	3, 522-531, 2001. 3. Li, Y. , Nemes, J.A. and Derdouri, A.A., Polym. Eng. Sci., vol. 41-8, 1399-1412, 2001. 4. Jones, A., 3 01 2 10 0 1 1 2 1 2 3 2 C S C F	41-

Identification of Mechanical Constitutive Equations


Nicolas Pirc

Florian Bugarin

Fabrice M Schmidt

Marcel Mongeau

BEM-based cooling channel optimization for 3D injection molding

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Today, around 30 % of manufactured plastic goods rely on injection molding, which is based on the injection of a fluid plastic material into a closed mold (Figure 1 [14] displays an illustrative injection molding process). The cooling time can represent more than 70 % of the injection cycle. Moreover, in order to avoid defects in the manufactured plastic parts, the temperature in the mold must be homogeneous. Thus, the design and the position of the cooling channels are crucial elements in the design of the mold. In order to decide the position and the shape of the cooling channels in the mold, designers commonly rely on experience and intuition within a costly trial-and-error design process. This manual design process becomes inadequate and unpractical for complex problems. This is particularly true nowadays with rapid prototyping processes such as layered design or selective laser sintering that enable manufacturers to build almost any desired shape of cooling channel geometry in the mold. As a consequence, designers need a more power-

INTRODUCTION

ful tool integrating the cooling analysis, its numerical simulation, and even optimization algorithms into the design process. We propose in this paper a practical methodology to optimize both the position and the shape of the cooling channels in 3D injection molding processes.

For the evaluation of the temperature, required both by the objective and the constraint functions, we must solve 3D heat-transfer problems via numerical simulation. Severe numerical methods such as Finite Element Method (FEM) [START_REF] Kolossov | 3D FEM simulation for temperature evolution in the selective laser sintering process[END_REF] or Boundary Element Method (BEM) [2] can be used for solving the heat-transfer problem. Silva [3] used 3D FEM software to model the injection molding cycle for complex geometries of molds. However, the computational burden of 3D FEM makes the integration of an optimization algorithm unpractical. Indeed, a FEM approach to the 3D heat-transfer problem imposes severe storage and CPU requirements, even if moderately complex industrial parts were to be targeted. An optimization algorithm requires numerous (and here expensive in CPU time) objective-function evaluations. Moreover, unless one knows a very good initial design, from one optimization iteration to the next the design is generally significally different. As a consequence, the capability of FEM to handle (small) shape perturbations for 3D mesh cannot be advantageously exploited here. BEM is a method that was popularized by Brebbia [6]. It is used in many applications, such as gas-assisted injection molding processes [4], or groundwater flow and mass transport problems [5]. BEM transforms domain integrals into boundary element integrals, involving therefore a discretization restricted only to the external and internal field boundary. An optimization method can therefore be envisaged to modify the position and also the shape parameters of the cooling channels in order to improve the cooling performance of the mold.

Park [START_REF] Park | Optimization method for steady conduction in special geometry using a BEM[END_REF] proposed such a BEM-based procedure for the position of the cooling channels relying on augmented-Lagrangian optimization method. However, his study is restricted to molds that are generalized cylinders, amenable to 2D molds. His optimization model minimizes the variations of the temperature distribution on the cavity surface with respect to the average temperature. As a consequence, his optimal configuration provides a uniform but high average temperature, leading a very long cooling time which is not desirable in the context of large-scale manufacturing. Mathey [8] also used BEM to solve the heat-transfer problems with a Sequential Quadratic Programming (SQP) [9] algorithm to improve mold injection cooling. She minimizes an objective function which is the weighted sum of two criteria. Her first criterion is the average temperature at the plastic-part surface. Her second criterion is the sum of the temperature variations with respect to the average temperature. Moreover, Mathey optimizes both the position and the shape of the cooling channels. However, her approach is also restricted to 2D molds (as for Park, BEM then reduces the dimension of the computation space from 2D to 1D).

Our contribution is threefold. First, we address 3D mold geometries with a BEM approach reducing the dimension of the computation space from 3D to 2D, avoiding full 3D remeshing: only the surface of the cooling channels needs to be re-meshed at each evaluation required by the optimization algorithm. Secondly, we propose a general optimization models that attempts at minimizing the desired overall low temperature of the plastic-part surface subject to constraints imposing homogeneity of the temperature. Thirdly, we demonstrate that our optimization methodology is viable with encouraging preliminary results on two semi-industrial plastic parts.

The paper is organized as follows. We detail in the next section the 3D heat-transfer problem that we must solve for every optimization evaluation of the temperature required by the optimization algorithm. Section 3 describes our overall optimization methodology. In order to validate our approach, we report encouraging preliminary results in Section 4. We conclude in Section 5.

THE HEAT TRANSFER PROBLEM

This section describes the heat-transfer problem that must be solved at every temperature evaluation required by the optimization algorithm: As shown in Figure 2 [10], after a few cycles the variation of the temperature in the injection mold in production can be considered to be quasi-stationary. Once the average temperature of the mold is stabilized, the cycle-averaged approach can predict well the overall performance of the cooling system. Thus, we can consider a stationary regime neglecting the transitory oscillations of the temperature. where ∆ is the Laplacian operator, T is a temperature vector. Multiplying this equation by a weighting function T*, and using Green's theorem, we obtain the well-known Somigliana's equation [6]. The following Green integral representation formula gives the value of T in terms of integral equations involving the fundamental solution T* of the Laplace equation:

C.T + Γ T.(∇T * .n).dΓ = Γ (∇T.n).T * .dΓ . (1)
Here, n is the unit normal at one element, Γ is the boundary of the domain, C is equal to 1 inside the domain Ω and to 0.5 on its boundary Γ. The weighting function T * and the flux q * are the fundamental solutions of the heat-problem, the so-called Green's functions [6]:

T * = -1 2.π.r and q * = -r.n 4.π.r 2 , ( 2
)
where r is the distance from the point of application of the concentrated unit source to any other point under consideration. Note that in (1) all the integrals are taken over the boundary of the domain. For the computation of the heat-transfer problem, the contour Γ is discretized and the integrals in the above equation are defined in terms of nodal values by means of interpolation functions. After reorganization of terms in (2) using the boundary conditions, we obtain a non-symmetrical linear system. We shall need to remesh the cooling channel surfaces using 2D elements, and to solve this system of equations at every optimization evaluation.

Figure 3 shows the boundary conditions on the mold.

λ.q = h(T a -T c ) ,
where λ is the conductivity, q is the heat flux, and T a is the ambient temperature. The following polymer properties are referenced in Table 1: ρ is the density and C p is the heat capacity. The flux density at the cavity surface, φ cavity , is computed from the cooling cycle time, t cooling , and the polymer properties (heat due to polymer crystallization is neglected) [START_REF] Jui-Ming | Multi-objective optimization scheme for quality control in injection molding[END_REF]:

φ cavity = Q t cooling .Γ cavity ,
where Q is the heat evacuated by the plastic part during one injection cycle. units polymer (PP) mold (steel 40cmd8s)

λ W m -1 K -1 0.63 34 ρ kg m -1 891 7800 C p J kg -1 K -1 2740 460
Table 1: Parameter values for heat-transfer problem Finally, remark that for our problem, we only need to mesh once the cavity surface and the external surface of the mold. Then, each time the optimization algorithm requires to evaluate the objective and the constraint functions, we must mesh the surface of the cooling channels. The output of the heat-transfer problem resolution is a set of temperature measurements, one for each surface element: cavity, external mold surface, and cooling channels. However, we only need the vector T of temperature measurements, {T i (x)} i∈S , at the cavity surface elements (S denotes the index set of the cavity surface elements).

OVERALL OPTIMIZATION METHODOLOGY

We first present in this section how we formulate our problem under a mathematical programming form. Then, we detail our overall computational methodology. In the sequel, x will denote the vector of optimization variables (position and shape parameters for the cooling channels). Since the output of the heat-transfer problem is function of x, we shall make explicit the dependence of the temperature measurements upon the position and shape parameters x : {T i (x)} i∈S .

Figure 4: The heat-transfer simulation coupled with the optimization algorithm.

Most practical optimization problems involve several (often contradictory) objective functions. It is the case here as one aims at minimizing the temperature of the plastic-part surface while minimizing the variation of the temperature along this surface. The simplest way to proceed in such a multicriterion context is to consider as objective function a weighted sum of the various criteria, which involves choosing appropriate weighting parameter values (more on multicriterion optimization can be found in [START_REF] Ehrgott | Multicriteria Optimization[END_REF]). An obvious alternative is to use one criterion as objective function while requiring in constraints maximal threshold levels for the remaining criteria. We choose here the latter approach because in our application context, it is more desirable to consider the temperature requirement as a constraint. Indeed, we do know a threshold level value for the maximal temperature variation under which any variation is equally acceptable. More precisely, we formulate our problem under the form:

min x T (x) subject to f (T (x)) ≤ 0 (3) g(x) ≤ 0, ( 4
)
where f is a real-valued function used to stipulate the uniformity-temperature constraint, and g(x) is a general vector-valued non-linear function. Remark that • stands for any norm such as the l 1 norm (sum of absolute values), the standard Euclidean (l 2 ) norm, or the max (l ∞ ) norm. The general constraints (4) represent any geometry-related or other industrial constraints, such as:

• upper/lower-bound constraints on the x i 's,

• keeping the cooling channels within the mold,

• technically-forbidden zones where we cannot position the cooling channels (for instance due to the presence of ejectors),

• constraints stipulating a minimal distance between every pair of cooling channels to avoid inter-channels collision.

Figure 4 shows the coupling between the thermal solver and the optimization algorithm.

APPLICATIONS

In this section, we present computational experiments on two semi-industrial plastic parts. For both applications, we specify the following implementation strategies. We use commercial software for the tasks represented by dashed-line boxes in Figure 4: • IDEAS for the initial 2D meshing of the cavity and the external surface (note that this meshing remains constant throughout the optimization process),

• Matlab's fmincon [9] Sequential Quadratic Programming subroutine for the optimization algorithm (with finite-difference computation of the gradients for these preliminary experiments-in future work, we intend to compute exact gradients).

We programmed every other task (represented by full-line boxes in Figure 4) in Matlab.

The system of linear equations involved in the heat-transfer problem is solved by the LApack [START_REF] Lapack | [END_REF] subroutine included in Matlab 7.0. We use here the l ∞ (max) norm for the objective function:

T (x) ∞ := max i∈S T i (x). ( 5
)
All numerical results reported in the sequel are obtained on a Macintosh 1.83 GHz Intel core 2 duo.

EuroTooling mold with straight cooling channels

For our first application, we consider a semi-industrial injection mold design for the European project: EuroTooling 21 [START_REF][END_REF]. The plastic part produced from this mold is a bended plate whose dimensions are shown in Figure 5. We want to optimize the position of the cooling channels. The channels are simple straight horizontal cylinders of constant length as illustrated in Figure 6. In order to simplify the presentation, in this application we only optimize the position of the cylinders (we could also consider optimizing the radius of the cross-section of the channels and the number of channels, etc.). As a consequence, three parameters suffice to describe the position of each cooling channel. For illustration purposes we consider here 8 cooling channels. The components of the optimization vector x are here the coordinates of each end point (X i , Y i ), i = 1 . . . 8. (in our application Z i is fixed). Thus, our problem involves 16 optimization variables. For this application, we use the constraint function F (T (x)) = T (x) -T (x) ∞ -σ, where the difference T (x) -T (x) stands here for the vector whose i th component is

T i (x) -T (x), T (x) := 1 |S| i∈S T i (x)
is the average of the T i 's, and σ is a user-defined temperature uniformity tolerance. We choose here σ = 4 • C.

In other words, we do not accept variations of the temperature that are above 4 • C with respect to the average temperature and this, everywhere on the plastic-part surface. The remaining constraints here simply stipulate lower/upper bounds and also forbid the grey zones displayed on Figure 6. The 2D meshing of the mold surfaces is displayed in Figure 10. The surface of each channel is discretized into 320 quadrangles. We use as starting point, a heuristic solution provided by an experienced engineer. On average, one objective-function evaluation (i.e. the computation of the temperatures) requires 86 seconds of CPU time. Since for these preliminary experiments we are content with computing gradients using finite-difference approximations, one optimization iteration involves 26 minutes of CPU time on average. Figure 7 displays initial (full lines) and optimal (dashed lines) positions of the cooling channels. Figure 8 displays the decrease of the objective-function value in terms of the number of optimization iterations. We observe on Figure 8 that the objective function is reduced by 90% of its initial value within the first four optimization iterations. Figure 9 shows the temperature distribution along the surface of the cavity mold before and after optimization (points A, B and B refer to positions on Figure 7). We observe that both the temperature variance and the average temperature decreased significantly. The temperature distribution after optimization can be seen on Figure 10

Plastic-cup mold with a helix cooling channel

We now report computational results on a 3D plastic part: a plastic cup. The height of the mold and that of the plastic cup are respectively 12 and 8 cm. The upper radius of the cup is 6 cm, its lower radius is 4 cm, and its thickness is 0.03 cm. In this preliminary study, we are content with considering a helix-shape cooling channel with circular cross section, with only three degrees of freedom: the section diameter of the cooling channel, the frequency of the turn and its overall diameter. The cooling channel is therefore defined in IR 3 using the single parameter algebraic equation:

X(t) = r s cos t Y (t) = r s sin t (6) Z(t) = t ,
where r s is the radius of the spiral. We discretize this curve into 63 straight-cylinder channels (Figure 12) with radius r c . Each of these cylinders is further discretized into 20 quadrangle elements as shown on Figure 13. The 2D meshing of the mold surfaces is displayed in Figure 11 The optimization variables here are:

• the geometrical parameters that control the shape of the helix, that is to say: the radius , r s of the spiral and the number, n, of helix turns,

• the radius r c of the circular cross section of the cylinder channels. Figure 14 shows the helix curve before and after optimization.

CONCLUSION

We proposed in this paper a practical methodology to optimize both the position and the shape of the cooling channels in industrial 3D molds. This was made possible through the use of the boundary elements method that avoids full 3D remeshing. Indeed, only the surface of the cooling channels needs to be remeshed to solve the heat-transfer problem involved each time the optimization algorithm needs to evaluate the temperature at the plastic part surface. Our optimization model can account for various ways to address the overall desired low-temperature criterion and the temperature-homogeneity constraint. Encouraging computational experiments on two semi-industrial plastic-parts, demonstrated the viability of our approach that is intended to be used as a decisionanalysis tool for designing new, original mold geometries. We are currently implementing exact gradients to replace finite-difference approximations in order to improve the efficiency of the optimization. Future work could allow more complex shape variations. Recall that we considered in this preliminary study only three degrees of freedom in order to test the viability of our optimization methodology. Other degrees of freedom for our plastic-cup application could for instance include the flux in the cooling channel and a vertical angle so that the helix runs along a cone (closer to the actual cup shape). Finally, further work could formulate our design problem as a topology optimization problem in order to offer even larger design freedom. The wide range of geometries we can tackle for the cooling channels will bring us to address more complex industrial molds. This will undoubtly yield challenging global optimization problems, as this can be anticipated based on the numerical experiments we reported.

Figure 1 :

 1 Figure 1: Temperature at the surface of the mold versus time

Figure 2 :

 2 Figure 2: Temperature at the surface of the mold versus time

Figure 3 :

 3 Figure 3: Boundary conditions

Figure 5 :

 5 Figure 5: Plastic part dimensions and one cooling channel

Figure 6 :

 6 Figure 6: Cooling channel parameters

Figure 7 :

 7 Figure 7: Initial and optimized (dashed lined) positions of the cooling channels

Figure 8 :

 8 Figure 8: Objective function versus iterations

Figure 9 :Figure 10 :

 910 Figure 9: Temperature measurements along the surface of the mold cavity before (black) and after (grey) optimization

Figure 11 :Figure 13 :

 1113 Figure 11: Meshing of the surfaces for the plastic-cup mold

Figure 14 :

 14 Figure 14: Helix curve before (black) and after (grey) optimization

Figure 15 :Figure 16 :

 1516 Figure 15: Objective-function value versus iterations

Acknowledgments

This study was conducted within the framework of the European project EuroTooling 21 [17].


Nadhir Lebaal

Fabrice Schmidt

Stephan Puissant

Optimizations Of Coat-Hanger Die, Using Constraint Optimization Algorithm And Taguchi Method

Keywords: Polymer extrusion, Response surface method, DoE, Optimization, Kriging Interpolation, Rem3D®

Polymer extrusion is one of the most important manufacturing methods used today. A flat die, is commonly used to extrude thin thermoplastics sheets. If the channel geometry in a flat die is not designed properly, the velocity at the die exit may be perturbed, which can affect the thickness across the width of the die. The ultimate goal of this work is to optimize the die channel geometry in a way that a uniform velocity distribution is obtained at the die exit. While optimizing the exit velocity distribution, we have coupled three-dimensional extrusion simulation software Rem3D®, with an automatic constraint optimization algorithm to control the maximum allowable pressure drop in the die; according to this constraint we can control the pressure in the die (decrease the pressure while minimizing the velocity dispersion across the die exit). For this purpose, we investigate the effect of the design variables in the objective and constraint function by using Taguchi method. In the second study we use the global response surface method with Kriging interpolation to optimize flat die geometry. Two optimization results are presented according to the imposed constraint on the pressure. The optimum is obtained with a very fast convergence (2 iterations). To respect the constraint while ensuring a homogeneous distribution of velocity, the results with a less severe constraint offers the best minimum.

INTRODUCTION

The design of dies for polymer extrusion often involves trial and error corrections of the die geometry to achieve uniform flow at the exit. Manual correction to die geometry is a time consuming and costly procedure. If the repartition channel in a flat die is not designed properly, the velocity at the exit of the flat die may not be uniform [START_REF] Wang | The flow distribution of molten polymers in slit dies and coat hanger die through three-dimensional finite element analysis[END_REF], and leads to a variation in the sheet thickness across the width of the die. Since a tight control on thickness is required for a high quality plastic sheet.

Often, the number of involved variables and their interactions prevent any optimization according to the trial and error corrections, because the number of evaluations needed may become very high. Hurez et al. [START_REF] Hurez | A New Design procedure for profile die[END_REF] used the trial and error corrections, to determine the optimal die shape which may allow to give a homogeneous exit flow distribution. Design of experiment, in particular the Taguchi method [START_REF] Montgomery | Design and analysis of experiments[END_REF], allows obtaining invaluable information on the important variables of the process in order to achieve the required goals. The effects of the various factors can be represented on graphs to support the discussion and to lead to identify the most sensitive to minimize the defects. Within this framework, we can mention Chen et al. [START_REF] Chen | Optimization of the Coathanger manifold via computer simulation and orthogonal array method[END_REF]. They showed, using the Taguchi method, that the operating conditions, the type of materials, and the geometry of the die have a great influence on the exit velocity distribution on the die.

Prior works in sheet die optimization have involved the use of lubrication approximations of the momentum equations [START_REF] Yu | A simple numerical approach for the optimal design of an extrusion die[END_REF][START_REF] Xiaorong | Optimal design for polymer sheeting dies[END_REF]. If the geometry is more complex, a flow channel can be approximated with simple geometric sections. Smith et al. [START_REF] Smith | Optimal design for polymer extrusion. Part I: Sensitivity analysis for nonlinear steady-state systems[END_REF][START_REF] Smith | Optimal design for polymer extrusion. Part II: Sensitivity analysis for weakly-coupled nonlinear steady-state systems[END_REF] modeled Newtonian and non-Newtonian isothermal flow in a coat hanger die using a generalized Hele-Shaw (HS) approximation, and optimized the die by minimizing pressure drop subject to exit flow uniformity being within a tolerance set. The sensitivities analysis needed for the Sequential Quadratic Programming (SQP) algorithm was calculated by direct differentiation and the adjoint method are compared, and simultaneous minimization of velocity dispersion subject to residence time variations is added using Broyden. Fletcher. Goldfarb Shanno (BFGS) algorithm and penalty function. The same author [START_REF] Smith | An optimisation-based approach to compute sheeting die designs for multiple operating conditions[END_REF] in order to optimize the geometry for two different materials at various temperatures, he use two optimization algorithms with constraint based on SQP and Sequential Linear Programming (SLP), and are compared. The optimization problem consists in minimizing the pressure loss in the die, with an imposed constraint so that a homogeneous velocity distribution is obtained on the outlet side of the die within an imposed tolerance. Network algorithms have been developed to optimize die designs [START_REF] Hurez | A New Design procedure for profile die[END_REF] but they are difficult to apply to arbitrary shapes. Michaeli et al [START_REF] Michaeli | Computer aided optimisation of extrusion dies[END_REF] have used a combination of finite element analysis for isothermal flow and flow analysis network to accelerate the iterative optimization process for the design of profile extrusion dies. To optimize the die geometry, they used respectively the evolution strategy algorithm and network theory. Sun et al. [START_REF] Sun | Optimization of a flat die geometry[END_REF] optimize a flat die using BFGS algorithm. A penalty function was introduced to enforce a limit on the maximum allowable pressure drop in the die.

The ultimate goal of this work is to optimize the die channel geometry in a way that a uniform velocity distribution is obtained at the die exit. For this end, we developed an automatic optimization algorithm based on a response surface method together with Kriging interpolation. We used the REM3D® software [START_REF] Silva | Viscoelastic Compressible Flow and Applications in 3D Injection Molding Simulation[END_REF] to compute a 3D flow in extrusion dies. This software takes into account strain rate and temperature dependence. The optimal design procedure is applied to a coat hanger die used to extrude thermoplastic sheets (Figure 1).

THE OPTIMIZATION BENCHMARK

In this paper, the geometry of a flat die is optimized for an Acrylonitrile Butadiene Styrene (ABS, Astalac EPC 10000). The extrusion simulation is carried out using the 3D computation software by finite elements REM3D®. The behaviors laws used in Rem3D® give an expression of the viscosity in function of the shear rate and temperature. The rheological parameters of the ABS are given in Table 1. Carreau Yasuda/WLF viscosity model is used to characterize the temperature and shear rate dependence [START_REF] Balasubrahman | Thermal control of melt flow in cylindrical geometries[END_REF]. It is written as: m-1 A flow of 50000 mm3/s was imposed on the entry with a temperature of 240°C, and the temperature of the die is constant and equals to 230°C.

Design Variables

The design of a flat die is based on the values of various geometrical parameters; the optimization variables (A, B and C) correspond to the geometry of the coat hanger die. The first variable A represents the depth of the channel repartition, the second variable B represents the opening of the channel repartition, the third one C represents the channel thickness (figure 1). During the process of optimization the variables A, B and C vary in a limited field. Each variable has its geometrical limitations which are respectively 30 < ,4 < 110 mm, 50 < B < 140 mm and 5 < C < 45 mm. During the optimization process, we developed a Matlab® based program which, allows us to treat the CAD and to change the geometry of the die automatically. In the first stage, in order to change the geometry we make a surfaced mesh. Since we know the nodes position, we seek to determine the coordinates of all the nodes (Ni), belonging to both surfaces SI and S2 (figure 2). For this, we use the diffuse approximation [START_REF] Lebaal | Rheological parameters identification using in-situ experimental data of a flat die extrusion[END_REF] with a linear interpolation. To change the geometry of the channel repartition we define the new surface (Sl k ), which corresponds to the new optimization parameters A k , B k , and C k . After that, we carry out a change of coordinates for all nodes which belong to the surface SI, so that they will be relocated on to the new Sl k surface which represents the new geometry. For S2 surface, a simple translation of the following nodes is applied there. The same procedure is applied on the other side of the die.

Objective And Constraint Functions

This optimization problem consists in determining an optimal geometry to homogenize the velocity distribution through the die exit, which corresponds to the minimum of the velocity dispersion (E). While preventing that the pressure "pressure loss" increases more than the pressure obtained by the initial geometry. We can also impose a more severe constraint on the pressure; this condition is translated by a constraint function (g).

min J(O)

Such that g •-P-(a*Po) (a*P 0 ) (2)

<0

where (J), the normalized objective function, the velocity dispersion (E), is defined as:

f E \_ N Z f\ v lY\ v (3) J)
where E 0 and P 0 are respectively the velocity dispersion and the pressure in the initial die, N is the total number of nodes at the die exit in the middle plane, v i is the velocity at an exit node, and v is the average exit velocity, is defined as:

1 N (4)
The constraint function (g) is selected in a way that it is negative if the pressure is lower than the imposed pressure, if not it will be positive. If a = 1; the pressure should not increase compared to the initial pressure. If a < 1; the pressure must be even lower compared to the initial pressure.

Study of the effects of the variables:

A Taguchi method [START_REF] Montgomery | Design and analysis of experiments[END_REF] is used to investigate the effect of the design variables (A, B and C). In the objective and constraint function we use the response obtained by a central composite design (CDD) [START_REF] Montgomery | Design and analysis of experiments[END_REF], which gives 5 levels per variable. The various levels are represented from 1 to 5 respectively giving the values of the minimum to the maximum of each variable. The effects of the various variables are represented on graphs to support the discussion and to lead to the identification of those influencing to minimize the defects.

Figure 3, illustrates the average effect of the design variables on the objective function at the various levels. From this figure we can note that the three geometrical variables A, B and C have considerable effects on the objective function, with strongly nonlinear variation according to the various levels. This is more particularly true for the variable A. We notice the variable "C" which represents the thickness of the repartition channel, has a more important effect at level 1 (the minimal thickness of the repartition channel). This effect represents a maximum value, corresponding to a worst exit velocity distribution. At the level 2, the effect represents the minimum value (or best exit velocity distribution). When C is at levels 2 to 5 the effect on the objective function is of same order as the other variables (A and B). It is noticed, from the same figure, that the variables A and B seem having means effects and oscillate compared to the various levels, which implies that a linear or quadratic mathematical model cannot represent effectively the objective function. As for their effects according to the various levels, we note that level 4 of variable A and level 5 of the variable B correspond at the minimum value of the objective function. Figure 4, illustrates the average effect of the optimization variables on the constraint function (pressure) at the various levels. It is also noted that the effect of the variable C on level 1 is most important. If the thickness of the repartition channel is minimal (CI) we need more pressure to extrude the plastic. For other levels of the variable C, we note that the effect is less important, and that the constraint function decreases gradually with the increase thickness of the repartition channel. Concerning the variable A, we note that the effect on the constraint function is average, and the minimum of the pressure corresponds to level 5. This level corresponds to the maximum depth of the repartition channel (figure 1). For this design, the extruded melt undergoes less contracting on the level of the repartition channel and flow out more easily, which can explain the pressure decrease. As for the variable B, we can see that its effect is small on the pressure (constraint), particularly for the levels 1, 2, 4, and 5. For level 3 of variables A and B, the value of the average effect increased because of the taking into account of the maximum value of the effect of the variable C on level 1.

f O 0.2 •0.2 -^-Variable C i i i i i i i --i - 1-1 r -r -r -i - i i i i i i i i i i i i i i i i i i i i i
^15 i 1 1 1 1 1 1 1 1 1 1-j|r-| 1 1 .5 ll l l l l l V "Q II I I I I I I 1 C ll l l l l l l I O ll l l l l l l 1 l4- 10 -i-^ 1 -i -i -i-i-^ ^ -i -V-i-i -£ II i i i i i i 1 •- II i i i i i i 1 **-_R i i i i i i i i i i i i i r ^ *ft m A1A2A3A4A5B1B2B3B4B5C1C2C3C4C5
Level of the variables of optimisation FIGURE 4. Effect of the variables on the constraint function.

Optimization procedure

The response surface method consists in the construction of an approximate expression of objective and constraint function starting from a limited number of evaluations of the real function. In order to obtain a good approximation, we used a Kriging interpolation described in next section. In this method, the approximation is computed by using the evaluation points by composite design of experiments.

Once the approximation of the objective function and the constraint function are built for each iteration. Since the successive evaluation of the approximate function does not take much computing time, we use the SQP algorithm to obtain the optimal approximate solution which respects the imposed nonlinear constraints. To avoid the fact of falling into a local optimum, and to respect the imposed nonlinear constraint, we use an automatic procedure which allows the launch of the SQP algorithm starting from each point of our experimental design. We take then the best approximate solution among those obtained by the various optimizations. Then, another approximate function is built by taking into account the weight function of Gaussian type which allows to slightly change the interpolation and makes the approximation more accurate locally (centred around the best minimum). Next, another minimization is carried out by the SQP algorithm with the initial point which represents the best optimum The iterative procedure stops when the successive points of the approximate function are superposed with a tolerance e=10" 6 . Finally, another evaluation is carried out to obtain the real response in the optimization iteration.

An adaptive strategy of the search space is applied to allow the location of the global optimum and, then, to readjust the research space by decreasing it by 1/3. A new design of experiment is automatically launched around the optimum. An enrichment of the interpolation is made by recovering responses already calculated, and which are located in the new research space. The iterative procedure stops when the successive points are superposed with a tolerance e=10 -3 .

Kriging interpolation

Kriging (27, 28), is an elegant and general interpolation method. Its origin goes up to Krige at the beginning of the 50 th , it was proved to be a very powerful interpolation technique. The Kriging, makes it possible to present the complexes function effectively (curved, surfaces...). This method is applied in our work to represent the response surface in an explicit form, according to the variables of optimization.

J(x) = p T (x)a + Z(x)

( The output responses from the function are given as:

F(x) = {f 1 (x\f 2 (x\-fA*)} (6)

From these outputs the unknown parameters a can be estimated:

a = (P 1 R-'PpP 1 R l F ( 7
)
Where P is a vector including the value of p(x) evaluated at each of the design variables and R is the correlation matrix, which is composed of the correlation function evaluated at each possible combination of the points of design:

R

^(*i, x i) •" ^(*i, x J i?(x n , Xl ) ... R(x n ,x n ) R 'j 4 (8) (9)
The second part in Eq. ( 5) is in fact an interpolation of the residuals of the regression model/? (x)a . Thus, all response data will be exactly predicted; is given as:

Z(x) = r T (x)P ( 10
)
Where r T (x)= {R(x,x 1 ),---,R(x,x n )}

The parameters /? are defined as fallow:

P = R l (F-Pa) (11)

RESULT AND DISCUSSION

A representative convergence history during an optimization process on the objective functions (J) is represented in figure 5. Two curves of results of the convergence of optimization are presented:

Results 1: A constraint on the pressure is imposed for which the pressure in the optimal die must be lower or equals 0.5*P0 {a = 0.5).

Results 2: The constraint imposed on the pressure is more severe. The pressure in the optimal geometry must be lower or equals 0.35*P0 (CC = 0.35).

Our algorithm of optimization clearly showed its capacity to obtain an optimal solution with a very fast convergence. The optimal solution is obtained in two iterations for both cases. The objective function is then reduced to: 95 % of its initial value, for the first result a = 0.5 . And 91% of its initial value, for the second result a = 0.5 We notice that the objective function for both results has strongly decreased in the first iteration, because of the precision of the Kriging interpolation. Nevertheless, to respect the constraint while ensuring a homogeneous distribution of velocity, we observe that result 1 offer the best minimum in the first and second iteration, because of the less severe constraint on the pressure compared to results 2.

We have represented the relative velocity dispersion obtained in the initial geometry and after optimization, for both results in figure 6. We see that the objective function obtained by imposing a severe constraint on the pressure (a = 0.35) is slightly higher of the objective function obtained with a less severe constraint (a = 0.5). The constraint on the pressure is represented on figure 7, in a form of (P/PO) for both results. We see that the constraint function is weaker with a more severe constraint (a = 0.35).

CONCLUSION

The applied method of optimization with the various strategies allowed us to obtain an optimal geometry in two iterations with a very satisfactory computing time for E.F calculation of the flow in three dimensions.

The optimal die gives a good exit velocity distribution while minimizing the pressure for the same flow. A qualitative comparison enables us to note that we were able to obtain with the optimal geometry a 95% profit of homogeneity compared to the initial geometry with a decrease of pressure of 52% (results obtained with a less severe constraint on the pressure). The use of the different constraints on the pressure allowed us to compare the various results obtained after optimization. The results are obtained at the second iteration. The algorithm of optimization clearly showed its robustness and its capacity to obtain an optimal solution with a very fast convergence. Total CPU time is of 3 days and 16 hours on a machine Pentium IV, 3 GHz, 1 Go RAM.

FIGURE 1 .

 1 FIGURE 1. Geometry of the flat die.

FIGURE 2 .

 2 FIGURE 2. Procedure for the change of the geometry.

FIGURE 3 .

 3 Effect of the variables on the objective function.

  with, p(x) = [p 1 (x), ..., p m (x)f, where m denotes the number of the basis function in regression model, a = fa v ..., a m ] is the coefficient vector the, x is the design variables, J(x) is the unknown objective or constraint interpolate function, and Z(x) is the random fluctuation. The term p (x)a in Eq. (5) indicates a global model of the design space, which is similar to the polynomial model in a Moving Least Squares (MLS) approximation. The second part in Eq. (5) is a correction of the global model. It is used to model the deviation from p (x) a so that the whole model interpolates response data from the function The construction of a Kriging model can be explained as follows:

FIGURE 6 .

 6 FIGURE 6. Velocity dispersion with initial and optimal die.

FIGURE 7 .

 7 FIGURE 7. Convergence history during an optimization run of the Constraint functions.

TABLE 1 .

 1 Rheological parameters.

	T]o Pa.s	m	CCo	Ts Pa
	716.997	0.15862	1	224063
	Al	A2[K]	Ts[K]	T ref [K]
	20.4	101.6	397.7	524.7

TABLE 2 .

 2 Summary of the results.

	Results of	Initial	Result 1	Result 2
	optimization	die	or = 0.5	or = 0.35
	"J"	1	0.049	0.09
	improvement	-	95%	91%
	"g" profit of pressure	1 -	0.48 52%	0.38 62%
	Variable A [mm]	36.5	96.6	103.08
	Variable B [mm]	112.6	76.15	72.74
	Variable C [mm]	25	25.47	31.4
	0		1	2
		Iterations
	FIGURE 5. Convergence history during an optimization
	process on the objective functions.

CP908, NUMIFORM '07, Materials Processing and Design: Modeling, Simulation and Applications edited by J. M. A. Cesar de Sa and A. D. Santos © 2007 American Institute of Physics 978-0-7354-0415-1/07/S23.00


M Bordival

F M Schmidt

Y Le Maoult

V Velay

Simulation of the Two Stages Stretch-Blow Molding Process: Infrared Heating and Blowing Modeling

Keywords: Stretch-blow molding process (SBM), heat transfer modeling, blowing simulation, G'sell law PACS: 44.05.+e, 44.40.+a, 83.60.St

In the Stretch-Blow Molding (SBM) process, the temperature distribution of the reheated perform affects drastically the blowing kinematic, the bottle thickness distribution, as well as the orientation induced by stretching. Consequently, mechanical and optical properties of the final bottle are closely related to heating conditions. In order to predict the 3D temperature distribution of a rotating preform, numerical software using control-volume method has been developed. Since PET behaves like a semi-transparent medium, the radiative flux absorption was computed using Beer Lambert law. In a second step, 2D axi-symmetric simulations of the SBM have been developed using the finite element package ABAQUS ® . Temperature profiles through the preform wall thickness and along its length were computed and applied as initial condition. Air pressure inside the preform was not considered as an input variable, but was automatically computed using a thermodynamic model. The heat transfer coefficient applied between the mold and the polymer was also measured. Finally, the G'sell law was used for modeling PET behavior. For both heating and blowing stage simulations, a good agreement has been observed with experimental measurements. This work is part of the European project "APT_PACK" (Advanced knowledge of Polymer deformation for Tomorrow's PACKaging).

INTRODUCTION

In a typical Stretch-Bow Molding (SBM) process, a Polyethylene Terephthalate (PET) preform is heated in an infrared (IR) oven to its forming temperature (around 100°C), and brought into contact with a mold of the desired shape. In such a process, the quality of the final bottle is closely related to heating conditions. Indeed, the preform temperature distribution has a strong effect on the blowing kinematic (stretching and inflation), and consequently on the thickness distribution of the final part. Temperature also affects the orientation induced by stretching, which, in turn, affects mechanical and optical properties of the bottle [1]. Temperature is therefore one of the most important parameters in SBM. However, its measurement remains a delicate task, especially in the thickness direction. Some experimental methods, such as IR thermography allows to measure the surface temperature during heating, but not its profile through the material thickness [2]. Recently the use of thermocouples inserted in the preform thickness was investigated in [3]. On the other hand, numerical methods are increasingly used. Researchers implemented models into commercial finite-element packages like ANSYS ® [4], FORGE3 ® [5], or developed their own software [6][7][START_REF] Martin | Proceedings of ANTEC[END_REF] with the aim of predicting the three-dimensional temperature distribution in the preform. Finally, some studies focused on the development of numerical optimization strategies for the SBM. Automatic preform shape optimization was proposed in [START_REF] Thibault | [END_REF], while optimization of heating system design was investigated in [START_REF] Bordival | 9th ESAFORM conference on material forming[END_REF]. The objective was to target a uniform temperature profile along the preform length.

The simulation of the blowing step has been also the subject of significant researches within the last two decades. Few studies focused on the feasibility of 3D temperature-displacement simulations [5,[START_REF] Wang | [END_REF]. But on the whole, researchers proposed 2D axi-symmetric models, with different material laws. A review is proposed in [5]. It can be noticed that the air pressure inside the preform is generally applied as a boundary condition, which can lead to unrealistic results [START_REF] Menary | 10th ESAFORM conference on material forming[END_REF]. Moreover, temperature distribution through the preform wall thickness is generally omitted.

In this work, a simulation of the two stage SBM process is proposed. The 3D temperature distribution of a rotating preform was computed taking into account all the process conditions, and the real oven design. In a second step, this temperature distribution (particularly through the wall thickness) was applied as initial condition for the simulation of the blowing step. For that, the finite element commercial package ABAQUS ® was used. Thanks to a thermodynamic model, the air pressure inside the preform is automatically calculated during simulation. Following sections focus on presenting each model.

PREFORM HEATING MODELING

In the SBM process, heating devices are often composed by a set of halogen lamps associated to aluminum reflectors. The preform translates through the oven, and is animated by a rotational movement to provide a uniform temperature along its circumference. Radiation emitted by the IR lamps is partially absorbed through the preform thickness, before being diffused in each space direction. Additionally, the preform tends to be cooled by air venting. In other words, preform reheating results from a combination between conductive, convective, and radiative heat transfers.

Heat Balance Equation

The evolution versus time of the preform temperature is governed by the following heat balance equation:

( ) r p q T k dt dT c ⋅ ∇ - ∇ ⋅ ∇ = ρ (1)
Where T = temperature, t = time, ρ = density, c p = specific heat, k = thermal conductivity, q r = radiative heat flux density. In order to solve this equation in 3D, a finite volume discretization is adopted. For that, the preform is meshed into hexahedral elements called control volumes. Equation ( 1) is integrated over each control volume and over the time, to obtain the following integro-differential formulation:

( ) ( ) dt d n q dt d n T k dt d t T c t r t t p Γ - Γ ∇ = Ω ∂ ∂ ∫ ∫ ∫ ∫ ∫ ∫ ∆ Γ ∆ Γ ∆ Ω . . ρ (2)
where Ω = control volume, Γ = surface of a control volume. Unknown temperatures are computed at the cell centre of each element. While the internal side of the preform is supposed to be adiabatic, the following boundary condition is applied to the external one:

( ) ( ) 4 4 ∞ ∞ - + - = ∂ ∂ - T T T T h n T k P PET P c P σ ε (3)
Where h c = natural heat transfer coefficient, ε PET = PET mean emissivity, σ = Stefan-Boltzman constant, T p = preform surface temperature at external side, T ∞ = ambient temperature. The method used for estimating PET mean emissivity is fully detailed in [2]. This boundary condition takes into account two types of thermal exchanges. The first one is due to the cooling by natural convection, the second one to the own emission of the preform. These exchanges are particularly important during the cooling stage.

Radiative Transfer Modeling

Over the spectral band corresponding to the IR lamps emission (0.38-10µm), PET behaves like a semi-transparent body. This involves that the radiative heat flux is absorbed inside the wall thickness of the preform, and can not be simply applied as a boundary condition. The radiation absorption must be taken into account through the divergence of the radiative heat flux, previously presented in the heat balance equation. This term represents the amount of radiative energy absorbed per volume unit; it is also more commonly called radiative source term. The computation of this source term can not be carried out without a precise understanding of radiative transfer properties, including its spectral and directional dependencies. Researchers proposed different numerical methods in order to compute the radiative source term, like raytracing [5] or zonal method [6]. The method used in this work is divided into two steps:

First of all, radiative heat fluxes reaching the preform surface are computed. For that, IR lamps are meshed into surface elements of which the contribution is taken into account via view factors computation. Moreover, IR lamps are assumed to behave like isothermal grey-bodies. Their emission is then defined by the Planck's law [START_REF] Modest | Radiative Heat Transfer[END_REF]. Finally, incident fluxes are calculated with the following equation:

( ) ( ) ( ) ti t i i ip T L S F q λ λ λ λ π ε ρ ∑ - = 1 0 (4)
Where ρ λ = PET reflexion coefficient, F ip = view factor between the lamp element i and the preform, S i = surface area of the lamp element, ε tλ = tungsten emissivity, L λ = Planck's intensity of the lamp i at the filament temperature T ti .

In a second time, the radiation absorption is computed according to the Beer-Lambert law (under the assumption of the non-scattering cold medium [START_REF] Modest | Radiative Heat Transfer[END_REF]): ( ) ( )

x q x q λ λ λ κ - = exp 0 ( 5
)
Where q λ (x) = spectral radiative heat flux density at the location x, q λ0 = incident spectral radiative heat flux density, κ λ = PET spectral absorption coefficient (in m -1 ).

Finally, the radiative source term is computed according to the following equation:

( ) λ κ λ κ λ λ λ d e q x q x r - ∆ ∫ - = ⋅ ∇ 0 (6)

Application -Results and Discussion

Software previously presented was used to simulate the reheating of a rotating preform with the processing conditions used on the laboratory blowing machine. The oven is composed of six halogen lamps (1 kW power), with ceramic and back aluminum reflectors. After 50 s heating, the preform is cooled down by natural convection during 10 s. The natural convection coefficient was calculated using the empirical correlation of Churchill and Chu [START_REF] Incropera | Fundamentals of Heat and Mass Transfer[END_REF]. Its value was estimated to 7 W.m -2 .K -1 . Percentages of nominal power of each lamp are reported TABLE 1. The preform rotating speed is equal to 1.2 rps. The preform used is 18.5 g weight, 2.58 mm thickness. The Material is PET TF9 grade (IV=0.74). An illustration is displayed FIGURE 1. Temperature measurements were performed in order to validate simulations. As it was demonstrated in [2], PET behaves like an opaque body over the 8-12 µm spectral band. For this reason, an AGEMA 880 LW IR camera, functioning within the long wave spectral band 8-12 µm, has been chosen. This choice makes possible to affirm that the camera measures a surface temperature. PET mean emissivity was also measured by following the protocol fully detailed in [2]. Its value is equal to 0.93. FIGURE 2 illustrates the external temperature distribution computed with the IR heating software, as well as the measured temperature cartography. In the aim of achieving more precise comparisons, the temperature profile along the preform length (at the end of the cooling step) is represented FIGURE 3. A good agreement between simulations and measurements can be observed, since the global error is less than 10%. FIGURE 4 illustrates the variation of temperature versus time on a single point, located at 47 mm from the neck of the preform (this point was chosen because it corresponds to the node located at the middle height of the mesh). This curve shows clearly the effect of the cooling stage. Indeed, it is interesting to notice that after 3 s of cooling (also called inversion time), temperature on internal side becomes higher than on the external one. This phenomenon can be easily

Pressure sensor A B

explained: while natural convection tends to cool the external side, the internal one is heated by heat conduction. In the SBM process, this point remains crucial. Indeed, there can be a significant difference between the inside and outside hoop stretch ratios. In order to ensure a good uniformity of the stress distribution through the thickness of the bottle, it is necessary to deliberately develop a non-uniform temperature profile throughout the preform before stretch and blowing. Finally, FIGURE 5 shows clearly that the temperature distribution through the thickness is not linear, but exponential. It can be seen that the temperature difference is around 4°C at the end of the thermal conditioning step. This value is of course strongly related to the cooling conditions. As a conclusion, results demonstrated the efficiency of the model developed in CROMeP. Infrared heating software remains a robust tool, allowing a better understanding of the effect of process parameters on temperature profiles, particularly through the preform wall thickness. It could also be used in order to optimize heating systems [START_REF] Bordival | 9th ESAFORM conference on material forming[END_REF]. However, a precise understanding of the effect of temperature on the blowing stage is necessary. For that, numerical model devoted to the simulation of the blowing stage was developed. Following section focuses on giving the key points about this model.

BLOW-MOLDING SIMULATION

Simulations of the SBM process were developed using the commercial finite element package ABAQUS ® . In this study, the objective is to simulate the process within the same conditions as on the CROMeP blowing machine, which means: simple mold for 50 cl water bottle and no stretch rod. A special attention was given to the measurement of each initial and boundary condition, namely temperature, air pressure, and heat transfer coefficient between the preform and the mold.

Boundary Conditions

As it was mentioned previously, the preform temperature distribution was measured and calculated in order to be applied as initial condition.

The heat transfer coefficient between the polymer and the mold was measured using a sensor developed for this study. Its peak value was estimated to 230 W.m -2 .K -1 , as illustrated FIGURE 6. The method used for this measurement is fully detailed in [START_REF] Bordival | 10th ESAFORM conference on material forming[END_REF]. This coefficient is of prime interest since it affects drastically the cooling time of the plastic bottle. Variation versus time of the air pressure inside the preform was measured using a Kulite sensor (FIGURE 1). As illustrated FIGURE 6, the air pressure follows typical variations. In the first time, the pressure increases sharply. As soon as the pressure is sufficient to blow the preform, air volume inside the bottle increases and consequently the pressure drops. While preform internal volume remains constant, the pressure reaches gradually its nominal value. This typical evolution of air pressure gives a good representation of the blowing kinematic. G. Menary [START_REF] Menary | 10th ESAFORM conference on material forming[END_REF] has shown that it is unrealistic to apply the pressure directly as a boundary condition. Indeed, the pressure drop would conduct to a deflation of the preform, and not to the rapid inflation observed experimentally. In this study, air pressure is not considered as an input variable, but is automatically computed thanks to the thermodynamic model "fluid element" available in ABAQUS ® . This model is based on the perfect gas law. Pressure measurements are only used for validating simulations.

Material Behavior

PET behavior was modeled with the following visco-plastic G'sell material model [START_REF] Gorlier | Caractérisation rhéologique et structurale d'un PET. Application au procédé de bi-étirage soufflage de bouteilles[END_REF]:

( ) ( ) ( )         - - =                     + - -       = t t a t m T T F with h T W F T k K 1 exp sinh exp exp 1 exp 2 1 0 ε β ε β ε σ & (7)
Where σ = equivalent Cauchy stress, ε& = equivalent strain rate, ε = cumulated strain, m = sensitivity to strain rate, (K,k 0 ) = consistence. This model takes into account both temperature and strain rate dependencies, as well as the strain hardening which appears for large deformations. It presents the advantage to be numerically stable and relatively easy to implement. However this phenomenological behavior law is reserved to a small range of temperature and strain rate. Moreover, it does not take into account the viscoelasticity of the material. Constitutive parameters have been identified using an inverse method (non-linear constrain algorithm called Sequential Quadratic Programming) from equi-biaxial tensile tests performed in Queen University of Belfast. The thermo-dependency was identified by [START_REF] Gorlier | Caractérisation rhéologique et structurale d'un PET. Application au procédé de bi-étirage soufflage de bouteilles[END_REF] from shear tests on PET T74F9. This model has been implemented within ABAQUS ® via a Fortran subroutine known as user creep.

Blow Molding FEM Model

In order to avoid long computation times, an axisymmetric model has been chosen. This approach is possible since both preform and mold designs are axisymmetric, as well as kinematic boundary conditions. The preform was meshed into 46 quadratic shell elements (96 nodes), with five integration points through its thickness in order to take into account the temperature gradient. The mold used is a prototype developed at CROMeP. It produces 50 cl bottle. This one has been assumed to be rigid and isothermal. Indeed, for one SBM cycle, its temperature increase is about 1°C [START_REF] Incropera | Fundamentals of Heat and Mass Transfer[END_REF]. In order to compute the heat transfer between the polymer and the mold, a coupled temperature-displacement model was chosen in ABAQUS ® Standard (implicit time integration scheme). The viscous dissipation was not calculated. However it could have an important effect on the preform temperature, and consequently on the blowing kinematic. As mentioned previously, no stretch rod is modeled. Finally, the contact between the preform and the mold is assumed to be stick. A good agreement is observed (around 15 % error on the mean thickness). We can notice that the measured thickness distribution is probably not optimal from an industrial point of view. This is due to the preform design used in this study, which is probably not adapted to this type of bottle shape.

Results and Discussion

Thanks to the thermodynamic model used in this study, it is possible to compare numerical and experimental blowing kinematics, by comparing the evolution of air pressure.

It can be seen FIGURE 9 that the pressure computed by the numerical model is not exactly the same as the measured one. However, tendencies are respected.

FUTURE WORK

Future work will aim to consolidate the model presented in this study. It is well known that PET behavior remains the key point for improving the model. Since the pressure curve gives a good representation of the blowing kinematic, it could be envisaged to couple the model to an optimization algorithm in order to identify automatically the constitutive parameters of the material law, by minimizing the difference between the measured pressure, and the calculated one. It is also crucial to investigate the influence of temperature distribution through the preform thickness on the blowing kinematic and on the thickness distribution of the final bottle. A sensitivity study can also be envisaged concerning the heat transfer coefficient mold/polymer, in order to prove its effect on the blowing.

FIGURE 1 .

 1 FIGURE 1. 18.5 g preform -PET T74F9 (IV=0.74).

FIGURE 2 .

 2 FIGURE 2. External temperature distribution after cooling -A: measured -B: simulated.

FIGURE 3 .

 3 FIGURE 3. External temperature profile along the preform length after 10 s cooling.

FIGURE 4 .

 4 FIGURE 4. Variation of temperature versus time.

FIGURE 5 .

 5 FIGURE 5. Temperature profiles through the preform wall thickness. Same location as FIGURE 4.

FIGURE 6 .

 6 FIGURE 6. Heat transfer coefficient and air pressure.

FIGURE 7

 7 FIGURE 7 illustrates the intermediate preform shapes versus time.

FIGURE 7 .

 7 FIGURE 7. Simulation of the preform shape evolution. Measurements of the thickness distribution of the final part were performed on bottles forming on the CROMeP blowing machine. Comparison with simulation results are illustrated FIGURE 8.

FIGURE 8 .

 8 FIGURE 8. Wall thickness distribution of the bottle.

FIGURE 9 .

 9 FIGURE 9. Computed and measured air pressure.

TABLE 1 .

 1 Process parameters of the IR oven

	P1	P2	P3	P4	P5	P6	t heat	t cool
	(%)	(%)	(%)	(%)	(%)	(%)	(s)	(s)
	100	100	18	5	50	100	50	10

ACKNOWLEDGMENTS

This study was conducted within the frame of 6th EEC framework. STREP project APT_pack; NMP -PRIORITY 3. www.apt-pack.com. Special thanks to Logoplaste Technology for manufacturing the preforms and Tergal Fibre for supplying the material, and QUB for giving tensile test results. Authors thank also V. Lucin for its contribution to this work.