Nadhir Lebaal
email: nadhir.lebaal@insic.fr

Fabrice M Schmidt
email: fabrice.schmidt@enstimac.fr

Stephan Puissant
email: stephan.puissant@insic.fr

Daniel Schläfli
email: daniel.schlaefli@maillefer.net

Optimization of coat-hanger melt distributors for the wire coating process

The objective of the rheological design of extrusion dies in the wire coating process is to distribute the melt around the conductor uniformly. If the manifold is not designed properly, the velocity at the die exit may be not uniform. However, in this paper we present an optimal design approach, this approach involves coupling an optimization routine with mesh generators and three-dimensional finite element simulation software. The objective is to obtain optimal coathanger melt distributor geometry to ensuring a homogeneous exit velocities distribution without increase the pressure at the die entrance. 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 and SQP algorithm in order to improve the exit velocities distribution by modify the manifold geometry using a simplified design model.

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. However manual correction to die geometry is a time consuming and costly procedure. CHEN et al. [1] showed using Taguchi method that, the operating conditions, the materials change and the die geometry have a great influence on the velocities distribution at the die exit. Wang et al. [2] investigate the effect of manifold angle and the contour of the manifold cross section on the flow distribution in the coat-hanger die using a threedimensional finite element code with isothermal flow of power low fluids. Winter and Fritz [3] presented a theory for the design of coat-hanger dies, with circular or rectangular section of the distribution manifolds. For the latter, provided a height aspect ratio of the manifold, the theory predicts material independence of the flow distribution. However, uniform distribution may not be achieved in practice. Considerable gains are to be expected from the use of adequate computational tools, for the numerical prediction of polymer melt flows through extrusion dies [4][5][6]. Nevertheless these computational tools are usually not user-friendly and require large CPU times. However, obtaining the optimum parameters is often the result of a tiresome work of trial and error corrections, during which, the various numerical solutions are tested and modified. The die complexity, the nonlinear relationship of the rheological behaviour complicate the design of die for polymer extrusion; further increases the difficulty of the die optimization problem and leading by this fact to a high number of numerical tests while the CPU time is often very high for three-dimensional simulations. Nevertheless, if the research of the quality criterion is made without specific methodology, it can easily lead to a solution which will not be optimal. The global automatic optimization of extrusion die design still remains a challenge and it is the purpose of this work to contribute to further developments of automatic die design tools. It was the need to have a design process less dependent on personal knowledge. This is also the main purpose of a die design code that is currently being developed, which main objective is to carry out the automatic search of the optimal flow channel geometry. The code consists of an optimization routine coupled with geometry and mesh generators and 3D computational software [7], based on the finite element method able to predict, the complex 3D flow and the temperature distribution. Because the quality criteria (objective and constraint functions) are implicit compared to the optimization parameters and their evaluations need a complete numerical analysis of the extrusion process, which requires a large computation times; the Response Surface Method [START_REF] Myers | Response surface methodology, Process and product optimization using designed experiments[END_REF][START_REF] Lebaal | [END_REF] (RSM) is adopted and coupled with an auto adaptive strategy of the research space to obtain the optimum parameters with low cost (acceptable computation times), and with a good accuracy.

2.Experimental data:

Die design:

The objective of this study is to find a design procedure for a distribution system that gives better uniformity of flow, i.e. uniform exit velocity and uniform average residence time over the die width. This uniformity should preferably be achieved with a die geometry that is independent of flow rate or polymer viscosity. Winter and Fritz [3] present a theory based on simplified flow model for the design of rectangular manifold to obtain a homogeneous exit velocity distribution. Given the polymer flow properties, the width of distribution system, 2b, the height of the slit h, and the volume flow rate

s v v bh Q 2 =
we can calculate the corresponding manifold contour lines y(x) as follows :

C x g x g x g x g Bb x y + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ + + - + - + = 1 ) ( 1 1 ) ( 1 ln 2 1 ) ( ) ( 1 2 3 ) ( (1)
with a function

( ) [ ] 1 2 1 ) ( ) ( - - = h x H x g (2)
The integration constant is equal to zero, C=0, given by the choice

h x H = ) ( at 0 = y
. The dimensionless parameter of the solution is defined as

b f h a B p = (3)
For a manifold of constant aspect ratio, a=W/H, the manifold has the side faces

( ) ( ) ( ) [ ] 3 1 h f a x b h x H p - = (4) ( ) [ ] 3 1 2 2 ) ( ) ( p f x b h a h x H a x W - = = (5)
Figure 1. Sketch of coat hanger distribution system with wide manifold and narrow slit flow region.

For a manifold of constant width H W >> , the die design is materials independent. The thickness variation of the manifold is then given by:

( ) ( ) W x b h x H - = (6)
The corresponding contour lines y(x) is given as follows :

( ) ( )

1 2 - - = W x b W
x y (7) In order to verify the design procedure, initial distributor is built according to the theory (eq.6-7) presented in [3] and tested. This coat hanger melt distributor has a diameter of 106 mm and a slit height of 3 mm, with a flat manifold of constant width W=50 mm. The initial manifold height is a little bit more than twice the slit height. Because of the high aspect ratio, the sidewall effects on the flow are small and almost independent of material changes. Therefore, the flat and wide manifold should result in a long (194 mm), but material insensitive distributor.

Materials:

A low density polyethylene (LDPE) referred Lupolen 1812D is selected for this example. A Carreau Yasuda/Arrhenius viscosity model is used to characterize the temperature and shear rate dependence of viscosity:

( ) a m a s 1 0 0 1 - ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + = τ γ η η η & (9) with ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = ref ref T T T T 1 1 exp ) ( ) ( 0 0 β η η (10)
The rheological properties (Tab. 1) are obtained by the data bases of Rem3D® (MatDB®). Table 1. Rheological parameters.

Material LDPE 0 η [Pa.s m ] m a s τ [Pa] T ref [K] β [K] 1812D
202278 0,35 1 10555 423 6156

Experimental measurements and numerical validation: 3.1 Experimental measurements:

The experimental measurements are established for LDPE Lupolen 1812D [10]. At the exit of the distributor, a flow separator divides the flow into 8 strips. The melt distribution is established by the measurement of the mass flow rate at some 36 seconds periods for each partial flow. It should be noted that a length of tools of 30 mm is used after the manifold length. Since the die and tools will have a more significant resistance to the flow, which will make it possible to improve the exit velocity distribution. By symmetry, the measurements of flow rate on sectors 1 and 5 remain unchanged. On the other hand, for sectors 2, 3, and 4 we respectively take the average of the measurements of flow established for the sectors ( Figure 2 illustrates the measured and calculated exit velocity distribution for the initial die design. The manifold of this initial die is of constant width and has an initial aspect ratio of h/w=0.11. According to the theory presented in [3], this design should be practically material-independent. However, distribution non-uniformity is important and perfectly predicted by the calculation (Figure 2). The velocity is more important in the sectors (1, 2, and 3) and weaker in the sections (4, and 5). However, the observed non uniform velocity distribution are not predicted by the theory; it is due to the already mentioned three-dimensional character of the manifold flow. Clearly, 3-D effects are important for such designs, but they are beyond the scope of the design theory. After these validation, the optimization procedure is introduced in order to improve the exit velocity distribution by modify the manifold geometry.

Optimization problem

Objective and constraint functions:

This optimization problem consists in determining an optimal geometry to homogenize the exit velocity distribution, which corresponds to the minimum of the objective function ( J ). While preventing that the pressure "pressure loss" increases more than the pressure obtained by the initial geometry; this condition is translated by a constraint function ( g ).

⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ ≤ - = = Φ 0 ) ( ) ( ) ( min 0 0 0 P P P g that Such E E J (15)
where E is the velocity deviation , is defined as follows:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = ∑ = N i i v v v N E 1 1 ( 16
)
where 0 P and 0 E are respectively the pressure and velocity deviation in the initial die, N is the total nodes number at the die exit in the mid plane, i v is the exit velocity, and v is the average exit velocity, is defined as:

∑ = = N i i v N v 1 1 (17)
The constraint function ( g ) is selected in a way to be negative if the pressure is lower than the pressure obtained by the initial die design (the pressure must be even lower compared to the initial pressure).

Optimization variables

For a diameter of 106 mm (b=166.5mm), and a slit height of 3 mm and initial manifold of constant width,

the contour lines y(x) and the thickness variation of the manifold are calculated starting from equations (6-7) given by Winter and Fritz [3]. During the optimization procedure, two variables witch represent the aspect ratio, a=W/H, and the shape factor, p f , will be optimized in order to ensure better exit velocity distribution. The new manifold contour lines y(x) and manifold thickness H(x) are calculated starting from equations (1-5).

Initially, the manifold width mm W 50 0 = is constant, and the shape factor p f is unity. During optimisation;

the aspect ratio a=W/H and the shape factor p f can vary respectively between:

10 5 . 0 ≤ ≤ k a and 1 0 ≤ ≤ k p f 4.

Sensitive analyses:

A Taguchi's design of experiments method [START_REF] Montgomery | Design and analysis of experiments[END_REF] is used to investigate the sensitive analysis of the design variables. A central composite design of experiment (CCD) [START_REF] Montgomery | Design and analysis of experiments[END_REF], is used to calculated the effect of each variable over all the field of research by studying the contribution of each variable on the responses variation "objective and constraint functions ". The CCD gives 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 graphical representations of factors effect at different levels are shown in Fig. 34. We note that the both variables (manifold aspect ratio "a" and shape factor " fp ") have more effect on the objective function and on the pressure. The optimum parameter level for minimum value of objective and constraint (pressure) functions are respectively a4 fp 4 and a2 fp (1-2).

In order to study the interactions between the two variables the graphical representations of factors interactions at different levels are shown in Fig. 5 (a andb). This figure represent the interactions between the two variables, by the average response with a , then fp in X-coordinate and by the value of their interaction, for the both functions (velocity distribution (Fig. 5-a) and pressure(Fig. 5-b)). We notes from this figure, that the effect on the objective and constraint function of the variable a is very strongly depends on the levels taken by the variable fp and reciprocally, witch implies that the interactions between the two variables is very strong. This strong interaction implies the non-linearity of the objective function.

Optimization procedure:

The RSM 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 [START_REF] Lebaal | NUMIFORM '07, AIP Conf. Proc[END_REF][START_REF] Kaymaz | [END_REF]. In this method, the approximation is computed by using the evaluation points by composite design of experiments. The Sequential Quadratic Programming b a (SQP) algorithm is used to obtain the optimal approximate solution which respects the imposed nonlinear constraints. After that, successive local approximations are built, in the vicinity of the optima by taking into account the weight function of Gaussian type. The iterative procedure stops when the successive optimum of the approximate function are superposed with a tolerance ε=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. During the progression of the procedure, the region of interest moves and zooms on each optimum. An enrichment of the interpolation is made by recovering responses already calculated, and which are located in the new search space. The iterative procedure stops when the successive points are superposed with a tolerance ε=10 -3 .

Optimization result and Discussion

A study of the effects and interaction of the optimization variables, show that the interaction between the optimization variables is more important. This, indicates that the non linearity of the function to be minimized is more important. The optimization example was carried out for LDPE 1812D and with a flow rate of 120 kg/h. By symmetry, only one half die is modelled. The optimization results are carried out on a machine Pentium D, 3.4 GHz, 3 Go RAM. 4. According to this table, we note an improvement of the objective function compared to the initial design of 93%. Furthermore, the constraint imposed on the pressure makes it possible to limit the increase in the pressure in the optimal die. We note that the constraint is respected and that the pressure decreased by 56 %. This is do the height manifold thickness.

The optimum solution is reached after two actualization of the research space. The two approximate objective functions are presented respectively in figure 67. The convergence history during the optimization run using SQP algorithm are shown for both approximations (1 st and 2 nd research space) respectively in figure 8 and figure 9. To avoid the fact of falling into a local optimum, and to respect the imposed nonlinear constraint, the best approximate solution among those obtained by the various optimizations is chosen. The true objective function is then reduced by 93% compared to their initial value. The calculated distributor performances are shown in fig 10. Differences between the two materials (LDPE) are predicted. However differences are minimal in term of velocity distribution. We noted in the seam figure that, after optimization, the dimensionless velocities distribution for the Lupolen 1812D are about unit for all die circumference (x). This implies a good velocity distribution compared by the initial die. Nevertheless, using the second material (LDPE 22H760) we observe week degradation of performance of a bout 5% in term of velocity distribution. We note in figure 12 that we obtain a good velocity distribution for two polymers. It is noticed that the constraint imposed on pressure is respected since the pressure did not exceed the initial pressure 17 MPa for the seam flow (120Kg/h). The pressure is higher in the case of Lupolene 1812D and it is weaker for the 22H760.

Because of the low aspect ratio, the sidewall effects on the flow are not insignificant and almost no independent of material changes. Therefore, the optimal flat and wide manifold result in a short (80 mm) compared to the initial manifold (193 mm), and the pressure drop decrease of about 50% compared to the pressure obtained by the initial die design (17.1 MPa).

Conclusions

The numerical algorithm presented in this work based on the global automatic optimization of extrusion die design shows promise as a tool capable of predicting optimal geometry and to have a design process less dependent on personal knowledge. The initial die obtained using the theory given by Winter and Fritz [3] gives a non uniform exit velocity distribution due the 3D effect. A comparison between the numerical results and experimental measurements

  where v and s i v are respectively the global average velocity and the average velocity for each sector, s j v is the velocity in each node j which belongs to the sectors i S and s N is the number of node for each sector on the die exit. The global relative deviation s g E is calculated by the sum of the relative deviations in each sector as follows:

Figure 2 .

 2 Figure 2. Calculated and measured velocity distribution.

Figure 3 .Figure 4 .

 34 Figure 3. Effect of factor levels on the objective functions.

Figure 5 .

 5 Figure 5. Interactions of factor levels on the objective (a) and constraint (b) functions.

Figure 6 .Figure 7 .Figure 8 .

 678 Figure 6. Approximate objective function (1 st research space)

Figure 9 .

 9 Figure 9. Convergence history of the objective function (2 nd research space)

Figure 10 .

 10 Figure 10. Exit velocity (dimensionless) distribution in the initial and optimal die.

Figure 11 .

 11 Figure 11. Pressures and exit velocity distribution in initial die design.

Figure 11

 11 Figure 11 presents the exit velocity distributions, and pressures in the initial geometry. The velocity distributions obtained by three-dimensional finite element model show the weakness of the design of one-dimensional model presented by Winter and Fritz [3]. The manifold flow can no longer be treated as a pure 1-D duct flow.

Figure 12 .

 12 Figure 12. Pressures and exit velocity distribution in optimal die design.

table 2 Table 2 .

 22 Process parameters

	2 and 8); (3 and 7); (4 and 6).
	3.2 Numerical simulations:
	The numerical simulations were carried out under
	Rem3D®. Due the symmetry, we model only half of
	the die cavity where the finite element model includes
	51512 nodes and 252593 elements. The process
	parameters ( v Q , Flow rate v Q [mm 3 /s]	Die wall temperature ( w T ) [°C]	Melt temperature (Tm ) [°C]
	47738	194	193
	For each sector, the dimensionless average velocity is
	defined by:

w T and Tm ), are referred to in

Table 3

 3 indicates the measured and calculated dimensionless average velocities using Rem3D® software for each sector at the die exit.

Table 3 .

 3 Results of Experimental measurements and numerical simulations:

	Tests	1
	Materials	LDPE Lupolène1812D
	Average velocities	a v exp	v	a cal
	S	1	1,098	1.101
	S	2	1,094	1.101
	S	3	1,071	1.059
	S	4	0,952	0.950
	S	5	0,785	0.789
	E	s g	52.6	52.23

Table .

 .

	4. A summary of the optimization results
	Optimization results	Initial die	Optimal die
	CPU time		18h40
	Iterations	0	2
	Objective function f	1	0.07
	Improvement of the velocity distribution	-	93%
	Constraint function P/P0	0.98	0.44
	Variable a	0.11	1.52
	Variable fp	1	0.36
	A summary of the optimization results is referred in
	table

Proceedings of the Polymer Processing Society 24th Annual Meeting ~ PPS-24 ~ June 15-19, 2008 Salerno (Italy)

indicates that the velocity distribution gives a perfect agreement. To improve the exit velocity distribution an optimization procedure is applied to predict the optimal manifold geometry by means of design model presented in [3]. The optimal 3D die extrusion is a short length compared to the initial die and gives a good exit velocity distribution while minimizing the pressure for the same flow.

Acknowledgements

The support of Maillefer SA is gratefully acknowledged.


B Cosson

F Schmidt

Y Le Maoult

NUMERICAL INFRARED HEATING OF SEMI TRANSPARENT THERMOPLASTIC USING RAY TRACING METHOD

Keywords: Heat transfer, Ray Tracing, Thermoforming

Thermoforming process include a heating stage of the thermoplastic specimen by infrared heaters. The knowledge of the temperature distribution on the surface and through the thickness of the preform is important to make good prediction of thickness and properties of the manufactured part. Currently in industry, the fitting of the process parameters is given by experience and it is expensive. Our objective is to provide tools that are able to simulate the heat transfers between the infrared heaters and the specimen in order to reduce the fitting cost and to control the qualities of the end product. We use the optical method called " ray tracing " (RTM) to simulate the radiative transfer. First, we compare the RTM with the shape factor method on a simple example: the heating of a square sheet by one infrared lamp. Then, we show the 3D heating stage of a rectangular specimen. In that case, we use 6 lamps to simulate an infrared oven. The RTM allows to compute a source term in the transient heat balance equation. We can use commercial codes to solve the heat balance equation.

INTRODUCTION

In this study we use two steps to solve the heat balance equation [START_REF] Cai | Dveloppement et applications de modèles d'échanges radiatifs par suivi de rayons[END_REF] of the heating stage of a Poly Ethylene Terephthalate (PET) plate by IR lamps. We compute the radiative source term Q (Eq.2) by using the ray tracing method [START_REF] Cai | Dveloppement et applications de modèles d'échanges radiatifs par suivi de rayons[END_REF], then we put the result in a finite elements software (here Comsol) as an input data. We choose ray tracing method because it is close to physics of IR heating. An other method, called radiosity, is often used in IR heating problem. This method is efficient for problem with opaque materials, but for problem with semi-transparent materials the radiosity method can not provide information inside the material.

ρC p dT dt = -∇ • (k∇ ¯T ) + Q (1)
where ρ, C p , k are respectively the density, the heat capacity and the thermal conductivity of PET, you can see F. Schmidt et al. for numerical values [START_REF] Bellet | Experimentai study and numericai simulation of the injection stretch/blow molding process[END_REF]. The radiative term Q, du to the IR lamps, is given by the equation ( 2).

Q = -K(T )M e -K(T )s (2) M = ∞ 0 M λ dλ (3)
K(T ) = ∞ 0 K λ M λ dλ ∞ 0 M λ dλ (4)
K λ is the spectral linear absorption coefficient of PET (Fig. 1), K(T ) is the integrated linear absorption coeffi- * Corresponding author: bcosson@enstimac.fr, 0033 563493085 cient that depends on the IR source (Fig. 2), M λ is the spectral exitance of the IR source given by the Plank's low.

VALIDATION OF RTM

RADIOSITY METHOD

In order to validate our software of ray tracing, we compare its results with an analytical solution (Eq.5), given by radiosity [START_REF] Person | Compilation of radiation shape factors for cylindrical assemblies[END_REF], in a representative case. We choose the configuration of a square sheet heated by one lamp (fig 4) of same length. The lamp is modelized by a cylinder. In both cases, we compute the power by unit of area (W.m -2 ) exchange between the lamp and the sheet.

F d1-2 = P S B -S 2Bπ            cos( Y 2 -B+1 A-1 ) + cos -1 ( C-B+1 C+B-1 ) -Y A+1 √ (A-1) 2 +4Y 2 cos -1 ( Y 2 -B+1 √ B(A-1) ) - √ C C+B+1 √ (C+B-1) 2 +4C cos -1 ( C-B+1 √ B(C+B-1 ) +H cos -1 ( 1 √ B )            (5) S = s/r ; X = x/r ; Y = y/r ; H = h/r (6) A = X 2 +Y 2 +S 2 ; B = X 2 +S 2 ; C = (H -Y ) 2 (7)
The geometrical variables (s, x, y, h) are defined on the figure 3.

Figure 3: Geometrical description of the cylinder-plate configuration

The equation 5 gives the illumination (W.m -2 ) of a differential element of plan illuminated by a cylinder. The result is given on figure 5 for the conditions described in figure 4. The power P of the lamp is 1000 W.

Lamp radiates in all the space. To use the ray tracing method, we have to discretize the radiation of the lamp and define ray properties : departure, direction and arrival. We use the Monte Carlo's method to define the ray departure and direction. We do not use regular discretization because the obtained result is discretization dependent. With the Monte Carlo's method we can control the error (Er F %, Eq.8) between the exact solution (Fig. 7.a) and thereby reduce calculating time. We can also control the convergence of the Monte Carlo's solution by looking the evolution of the mean of the illumination of the plate (Fig. 7

.b).

Er F % = 100

F radisity d1-2 -F raytracing d1-2 2 dΩ F radisity d1-2 2 dΩ (8)
For 1.2 million of rays, the error between the solutions given by radiosity and ray tracing method is smaller than 5%. The calculating time is about 90s with an Intel Centrino CPU (2.6 GHz). The ray tracing method is validated in the case of a plate illuminated by one lamp. In the next section, we compare numerical results to experimental results and we propose a technique to reduce the calculating time in the case of 3D computation.

COMPARISON WITH EXPERIMENT

In this section we compare the results given by our ray tracing software, coupled with Comsol (FEM), to experimental results given in [START_REF] Champin | Modélisation 3D du chauffage par rayonnement infrarouge et de l'étirage soufflage de corps creux en PET[END_REF] for the heating stage of a rectangular plate by an IR lamp (Fig. 8). In a first step, we compute the radiative source term (Eq.2), by using the ray tracing method, in all the plate (3D). Then, we use the result as an input data in a FEM software (here Comsol) to solve the heat balance equation (Eq.1). The figure 9 shows the source term for y= 0 in W.m -3 . The radiation of the IR lamp is absobed in few µm in the plate. For the chosen example, the heating stage of the plate is in two steps: in first the face of the plate that is in front of the lamp is heated, during 65 s, then the lamp is turning off and the temperature, in the thickness, is homogenised by conduction, during 25 s. The figure 10 shows the evolution of the temperature in two points of the plate: one in the middle of the face that is in the front of the lamp (front face) and the second is at the opposite (rear face). On the same figure, there are the numerical and the experimental results. At the end of the heating stage (90 s) the difference between simulation and experimentation is under 5 %. The numerical results are the dotted line, the blue one is for the front face and the green one is for the rear face. The grey lines are for the experimental results. The boundary conditions used for the simulation are given in the table 1. The computation time for this simulation is about 15 mn (Ray Tracing) for the source term and 3 mn (Comsol) for the temperature with a 2.6 Ghz CPU. The computation time of the source term is proportional to the number of lamps. In industrial oven you can have over 6 lamps and for this configuration the computation time is up to 1.5 hour. This time is acceptable for one simulation, but it can

TIME OPTIMISATION

In the simulation of the heating stage of a plate (or preform in a most general case), the calculation of the radiative source term is very time costly. In order to reduce the computing time cost, we propose to use interpolation for the calculation of the radiative source term. To realize this interpolation, we do a fixed number of direct calculation of the source term (Fig. 11). In our example, we can We use the interp2 function in Matlab, with a cubic interpolation, to realize the interpolation of the radiative source term. In order to validate this technique, we compare the result of a direct calculation and the result given by interpolation for arbitrary selected positions, in all our tests the error (Er Q %, Eq.9) between both results (Q interp and Q direct ) is under 3 %. So, with this technique the calculation time of the radiative source term comes from 15 mn to 3 s.

Er Q % = 100 (Q interp -Q direct ) 2 dΩ (Q direct ) 2 dΩ (9)

CONCLUSIONS

In this study we use the ray tracing method to compute the radiative source term used in the heat balance equation. The ray tracing method is validated with two different ways : comparison with analytical solution and comparison with experimental data. In both cases the error is under 5 %. The error can be controlled by using the Monte Carlo's method.

We propose a simple technique to reduce the computation time cost of the ray tracing method in the case where we need a lot of iterations, like in optimisation of forming process. By using this technique, we introduce an error smaller than 3% and reduce the computation time from 15 mn to 3 s. Next steps of this study are : the consideration of reflectors, the application to the heating stage of preform for the stretch blow moulding process.

Figure 1 :

 1 Figure1: Spectral linear absorption coefficient K λ[START_REF] Champin | Modélisation 3D du chauffage par rayonnement infrarouge et de l'étirage soufflage de corps creux en PET[END_REF]

1 )Figure 2 :

 12 Figure 2: Integrated linear absorption coefficient K(T )

Figure 4 :Figure 5 :

 45 Figure 4: Square plate with 1 lamp

Figure 6 :

 6 Figure 6: Illumination of the square plate (W.m -2 , ray tracing)

Figure 7 :

 7 Figure 7: (a) Error between radiosity and ray tracing methods, (b) Convergence of the ray tracing methodbe unacceptable if you want to do an optimisation of the heating stage. In the next section we propose a technique to reduce the computation time for a lot of simulations.

Figure 8 3 Figure 9 : 3 )Figure 10 :

 839310 Figure 8: Experimental configuration

Figure 11 :

 11 Figure 11: Position of the lamps in the plan x,z

Table 1 :

 1 Simulated boundary conditions

ACKNOWLEDGEMENT

This research was supported by Toshiba Lighting France.


M Bordival
email: bordival@mines-albi.fr

F M Schmidt

Y Le Maoult

B Cosson
email: cosson@ensm-douai.fr

A RAY TRACING METHOD TO SIMULATE THE INFRARED HEATING OF SEMI-TRANSPARENT THERMOPLASTICS

Keywords: Ray Tracing, Heat transfer modelling, Injection stretch blow moulding, Infrared heating, PET

This paper focuses on describing RAYHEAT, a software based on the ray tracing method, developed in order to simulate the InfraRed (IR) heating of semi-transparent polymers. In this study, RAYHEAT is used to simulate the IR heating step of a PET preform for the Injection Stretch-Blow Moulding (ISBM) process. The general principle of the method is to discretize, into a set of rays, the radiative heat flux emitted by halogen lamps, then to follow these rays inside the oven while they are not fully absorbed. The ray tracer computes the optical path of each ray -accounting for specular or diffuse reflections, refractions, etc… -from its emission point, and throughout the preform thickness. PET is assumed to behave like a non-scattering cold medium. Thus, the radiative heat flux absorption inside the preform is computed according to the Beer-Lambert law. Finally, the distribution of the radiative source term is calculated in the preform. In a second step, the radiative source term is applied as an input data in the commercial finite element software ABAQUS ® , in order to calculate the 3D temperature distribution in the preform. The source term is assumed to be time dependent in order to account for the preform movement throughout the IR oven. This method provides relatively small computation times, while keeping the memory requirements down to a minimum. Numerical results have been compared with temperature measurements performed on an in-lab IR oven. The model simulates suitably the infrared heating stage, and provides accurate predictions of the temperature distribution in the preform. The relative error between the temperature calculated by RAYHEAT, and the measured temperature, is less than 5%.

INTRODUCTION

A large part of bottles intended to beverage market is manufactured using one-step Injection Stretch-Blow Moulding (ISBM) process. ISBM involves an InfraRed (IR) heating stage to condition the PET preform to the appropriate temperature distribution (around 100°C). Then, the preform is stretched-blown inside a mould to get the final bottle. The heating conditions, which control the preform temperature distribution, strongly affect the final properties of the bottle (mechanical, optical and barrier properties). Within the last ten years, significant researches have been carried out on the simulation of the IR heating stage. A recent literature review on this subject has been presented by Bordival et al. [START_REF] Bordival | Optimization of Preform Temperature Distribution for the Stretch-Blow Molding of PET Bottles: Infrared Heating and Blowing Modeling[END_REF]. PET behaves like a semi-transparent body over the IR spectral range, resulting in major challenges with regard to radiative heat transfers modelling. Different approaches are presented in the literature in order to compute the radiative absorption through the preform thickness. The most classical is the view factors method [START_REF] Monteix | Experimental study and numerical simulation of preform sheet exposed to infrared radiative heating[END_REF][START_REF] Huang | Temperature profiles within reheated preform in stretch blow molding[END_REF][START_REF] Yousefi | Simulation of the mobile preform reheat in injection stretch blow molding process[END_REF]. This method consists in calculating the radiation on the preform outside surface, and then to assume that the radiative heat flux is unidirectional through the preform thickness. More recent works investigated the efficiency of the zone method [START_REF] Michaeli | FE-Analysis of the Two-Step Stretch Blow Molding Process[END_REF] or the ray tracing method [START_REF] Champin | 3D Finite Element Modeling of the Blow Molding Process[END_REF]. The last one presents the advantage to account for a while variety of optical effects, such as specular reflection or refraction. In addition, ray tracing enables to take into account most of constitutive elements of an IR oven like multiple lamps (various geometries) and reflectors (ceramic or metallic made). In this work, we propose a numerical modelling of the IR heating stage. The temperature distribution of the preform is calculated using a two-steps simulation. First, in-lab software, called RAYHEAT, is used in order to compute the absorption of the radiative heat flux inside the preform. Then, the results are applied as an input data in the finite elements commercial package ABAQUS ® . The method is applied in order simulate the IR heating stage of a rotating preform using the in-lab blowing machine set-up. A numerical validation is performed using temperature measurements. Following sections focus on presenting the governing equations of the model.

HEAT TRANSFER MODELLING

HEAT BALANCE EQUATION

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

( ) r p q T k t T c ⋅ ∇ - ∇ ⋅ ∇ = ∂ ∂ ρ (1)
where T = temperature, t = time, ρ = density, c p = specific heat, k = thermal conductivity, q r = radiative heat flux density. The inside surface of the preform is assumed to be adiabatic, while the following boundary condition is applied to the outside surface:

( ) ( ) 4 4 ∞ ∞ - + - = ∂ ∂ - T T T T h n T k P PET P c P σ ε (2)
where h c = convection heat transfer coefficient, ε PET = PET's mean emissivity, σ = Stefan-Boltzman constant, T P = outside surface temperature, T ∞ = air temperature inside the oven. The boundary condition given by Equation (2) accounts for two types of heat transfer. The first one is due to the cooling by natural convection, the second one to the preform own emission. These heat losses have a critical effect, especially throughout the cooling stage.

RADIATIVE SOURCE TERM

Over the spectral range corresponding to the emission of IR lamps (0.35-8 µm), PET behaves like a semitransparent body. This involves that the radiative heat flux is absorbed inside the wall thickness of the preform, and cannot be simply applied as a boundary condition. The radiation absorption must be taken into account through the divergence of the radiative heat flux (Equation 1), also called "radiative source term". The radiative heat flux is required to compute the radiative source term. Per definition, the radiative heat flux density is given by:

( ) ( ) ∫ ∫ ∞ Ω Ω Ω = 0 4 , π λ λ d d s I s q r (3)
where I λ = spectral intensity at the point s, along the direction Ω. The variation of the spectral intensity is governed by the Radiative Transfer Equation (RTC), given for a non-scattering medium by [START_REF] Modest | Radiative heat Transfer[END_REF]:

[ ] ) , ( )) ( ( ) , ( Ω - = Ω s I s T I ds s dI λ λ λ λ κ o (4)
where I°λ = Planck's spectral intensity at the preform temperature T, ds = optical path, κ λ = PET spectral absorption coefficient. Equation (4) has got an exact solution, given by [START_REF] Viskanta | Heat Transfer in Semitransparent Solids[END_REF]: In Equation ( 5), the first term of the right member represents the fraction of the incident intensity (emitted by the IR oven) which is absorbed by PET. The second term represents the fraction of the intensity, emitted by PET, which is absorbed by PET (due to the own emission). Throughout the heating stage, the preform temperature (approximately 400 K at the end of heating) is very low compared to the IR heater temperature (typically around 2400 K). In this case, the second term of Equation ( 5) can be neglected. This assumption leads to:

( ) ( ) ( ) ( ) ( )
( ) ( ) s e I s I λ κ λ λ - Ω = Ω , 0 , (6)
Equation ( 6) is called Beer-Lambert law [START_REF] Modest | Radiative heat Transfer[END_REF]. Then, by assuming that lamps emission is Lambertian, the radiative source term is given by [START_REF] Erchiqui | Modlisation par éléments finis du chauffage infrarouge des membranes thermoplastiques semi-transparentes[END_REF]: [START_REF] Modest | Radiative heat Transfer[END_REF] where M λ (0) = incident spectral emissive power. The radiative source term distribution is assumed to be temperature-independent. Therefore, it can be computed and then applied as an input data in the finite element software ABAQUS ® in order to solve the heat balance equation.

( ) ∫ ∞ - - = ⋅ ∇ 0 . ) 0 ( λ κ λ κ λ λ d e M s q s r

RAY TRACING METHOD

Ray tracing is widely acknowledged as one of the most accurate method used to simulate radiative heat transfers in semi-transparent media. The general principle of the method is to discretize, into a set of rays, the radiative heat flux emitted by halogen lamps. Then, each ray is followed inside the oven, from its emission point and throughout the preform thickness. The ray tracer computes the optical path of each ray, accounting for specular or diffuse reflections, refractions, etc… In our ray tracing software, assumptions are made for the different optical properties of lamps, reflectors and preforms. Those assumptions are referenced in table 1 for the following properties: emission, absorption, reflection and refraction. Reflection and refraction are averaged in order to reduce the number of ray stored in each calculation. In fact, for one ray that comes from the lamp and contained all the spectral information, if spectral reflection (or refraction) is computed, one ray by spectral band (an infinity for an exact model) has to be created for each air-PET interface crossing. Only the lamp filament is taken into account. Filaments are modelled by equivalent cylinders, the spiral form is neglected. Tungsten filaments are assumed to be Lambertian grey bodies. This assumption provides the definition of ray direction vectors (Figure 1) for rays coming from the filament.

Figure 1: Ray definition for ray tracing

The direction vector is defined by two parameters: θ and φ, respectively defined in the ranges [0, π/2] and [0, 2π]. Notations are illustrated on Figure 2. The computation of θ and φ has a strong effect on the ray tracing accuracy. Determinist discretization of the emission space could lead to errors due to the ray effect [START_REF] Li | Ray effect in ray tracing method for radiative heat transfer[END_REF]. To avoid that, we have chosen to compute θ and φ according to stochastic variables:

2 1 2 ; ) ( arcsin R R π ϕ θ = = (8)
where R 1 and R 2 are independent uniform stochastic variables in the range [0; 1]. The direction change of a ray that cross a PET-air interface is given by the Snell-Descartes law: where n air (= 1) is the refractive index of air. Notations are illustrated on Figure 2. In order to validate numerical simulations, IR heating trials have been performed on an in-lab IR oven. The oven is composed of five halogen lamps (2 kW nominal power), with ceramic and back aluminium reflectors (see Figure 3). Percentages of nominal power of each lamp are reported in Table 2. After 25 s heating, the preform is cooled down by natural convection during 10 s. No ventilation system is used. The natural convection coefficient was calculated using the empirical correlation of Churchill and Chu [START_REF] Churchill | Correlating Equations for Laminar and Turbulent Free Convection from a Vertical Plate[END_REF]. Its value was estimated to be 7 .5 W.m -2 .K -1 . The preform used is 27 g weight, 3.85 mm thickness. Its rotational speed is equal to 1.1 lap/s. The grade of polymer is PET INVISTA 1101 (IV = 0.8).

Figure 3: In-lab infrared oven

Temperature measurements have been performed using an AGEMA 880 LW IR camera, functioning within the long wave spectral range 8-12 µm. Over this spectral band, PET behaves like an opaque body [START_REF] Monteix | Quantitative infrared thermography applied to blow moulding process : measurement of a heat transfer coefficient[END_REF]. Thus, the camera measures a surface temperature.

PET THERMAL/RADIATIVE PROPERTIES

PET radiative properties were measured according to the protocols defined by Monteix et al. [START_REF] Monteix | Quantitative infrared thermography applied to blow moulding process : measurement of a heat transfer coefficient[END_REF]. Measurements were performed on PET INVISTA 1101 samples using a Perkin Elmer FTIR spectrometer over the range 2.5-25 mm. PET thermal properties are assumed to be temperature-dependent, which is particularly true for the heat capacity that sharply increases above the glass transition temperature.

NUMERICAL PARAMETERS

The preform is meshed into 5220 rectangular linear elements (6000 nodes), with ten nodes in the thickness direction. For each lamp, 1.5 million rays are followed, for a total of 7.5 millions rays for the entire oven. The computational time of the source term is equal to 90 min CPU (T9500 2.6 GHz 3Go RAM). Temperature calculation was performed using ABAQUS ® standard 6.7, with a fixed time increment equal to 0.25 s. The computation time is equal to 12 min CPU.

RESULTS AND DISCUSSION

Figure 4 illustrates the radiative source term distribution versus time. We can observe that the source term has been defined as a function of time in order to account for the preform rotation. Figure 5 illustrates the external temperature distribution along the preform height, at the end of the cooling step. We can observe a good agreement between simulations and measurements, since the mean relative error is equal to 1.7 %. Figure 6 illustrates the variation of external temperature versus time on a single point, located at 40 mm from the neck of the preform (middle-height). The agreement is fair, especially during the cooling stage, which indicates the convection heat transfer coefficient is well estimated.

CONCLUSIONS

We have proposed in this work a modelling of the IR heating stage for the ISBM process. The ray tracing software RAYHEAT provides accurate predictions of the 3D radiative source term distribution in the preform. The coupling between RAYHEAT and ABAQUS ® allows the computation of the preform temperature, accounting for the preform movement throughout the oven. A numerical validation has shown the efficiency of the method, since the relative error between simulations and experiments is less than 2%. Future work will aim to further improve the ray tracer in order to decrease the computation times.

Figure 2 :

 2 Figure 2: Ray refraction at air-PET interface crossing

Figure 4 :Figure 5 :

 45 Figure 4: Source term distribution versus time

Figure 6 :

 6 Figure 6: External temperature versus time (mid-height)

Table 1 :

 1 Assumptions on optical properties

		Emission Absorption Reflection	Refraction
	Lamps	Spectral isotropic	none	none	none
	PET	Averaged	Spectral	Specular and averaged	averaged
	Ceramic reflector	none	opaque	diffuse and averaged	none
	Metallic reflector	none	opaque	Specular and averaged	None

Table 2 :

 2 Lamps setting for IR heating trials

		L 1	L 2	L 3	L 4	L 5
	Power (%)	100 100 20	5	60

ACKNOWLEDGEMENT

Special thanks to SIDEL company (http://www.sidel.com) for manufacturing the preforms and for their collaboration. Authors also thank Mr Plantamura for his contribution to this work.