S Brosillon 1 '

M-H Maner0

J N Foussard

ADSORPTION OF ACETONE/ HEPT ANE GASEOUS MIXTURES ON ZEOLITE CO-ADSORPTION EQUILIBRIA AND SELECTIVITIES

Keywords: VOC removal, zeolite, fixed bed, co-adsorption, breakthrough curves

In the present study, co-adsorption of gaseous n-heptane and acetone is studied. The adsorbent is a commercial hydrophobie zeolite. Breakthrough curves and isotherms are measured for three different acetone/heptane mixture ratios and for pure components. Zeolites showed a greater adsorption capacity for acetone alone than for heptane alone. For the mixtures, the majority component was the most adsorbed. Zeolites have a good selectivity for acetone even though acetone is more volatile than heptane. Specific adsorption sites are reserved for acetone adsorption and are not accessible to heptane. The polarity of a component leads to electrostatic interactions with cation exchange. On non-specific adsorption sites, acetone can be displaced by heptane, which is Jess volatile. Moreover, experiments show that selectivity depends on the ratio of the mixture. For co-adsorption of a mixture of polar and non-polar components, selectivity depends on relative polarity, mixture ratio, and boiling point.

INTRODUCTION

Since new regulations on volatile organic compound (VOC) emissions have been decided, industrial cornpanies have become more and more interested in ways to reduce their pollution. Conventional engineering processes can be used to remove VOCs. These removal processes can be classified into two categories: destructive techniques and recovery techniques [START_REF] Manero | Industrial air pollution : removal of dilute gaseous vapor[END_REF]. ,.---0 (heptane or acetone), the breakthrough curves are dispersed;

Destructive techniques

Combustion is the most commonly

(ii) For the minority compound (heptane or acetone), the breakthrough curves are superimposed or very close. These tendencies are connected with the isotherm shape. The initial concentrations for the rich mixture are probably not on the plateau of isotherm, whereas those for the poor mixture probably are. This is why the breakthrough curves are superimposed.

Single component

The adsorption isotherms at 298 K for heptane alone and acetone alone are shown in Figure 5. Even if acetone is more volatile than heptane, HISIV 3000 adsorption capacities plateau for a heptane-rich mixture, equimolar mixture, and heptane-poor mixture represent 94%, 32% and 31% respectively of the maximal amount adsorbed in the case of heptane alone. Adsorption of a heptane-rich mixture and pure heptane adsorption have similar results. Equimolar conditions seem not to be favourable for heptane as the amount adsorbed is close to the amount adsorbed for heptane-poor mixture. The shapes of the isotherms can be explained by breakthrough results. On the one hand, dispersed breakthrough curves lead to a round isotherm; whereas superimposed breakthrough curves involve a flat isotherm.

For acetone (Figure Sb), the amount adsorbed at the plateau for an acetone-rich mixture, an equimolar mixture, and a heptane-poor mixture represent respectively 45%, 23% and 22 % of the maximal amount adsorbed in pure acetone adsorption. In this case the acetone isotherm for an acetone rich mixture is different from pure acetone isotherm. Adsorption of an acetone-rich mixture is not similar to adsorption of acetone alone. In the same way as for heptane, the adsorption isotherm of the equimolar mixture and the acetone-poor mixture are very close.

In Figure 6, these results are plotted on three different diagrams, each representing a given pair: heptane-rich, acetone-rich and equimolar mixture. As seen before, f o r a heptane-rich mixture, the amount of heptane adsorbed is close to the amount adsorbed when only heptane is present. Moreover it can be seen from Figure 6 that the amount of acetone adsorbed is not equal to zero. So there are sites accessible only to acetone where heptane is unable to desorb it. For an acetone-rich mixture (Figure 6b), as seen before, acetone is better adsorbed than heptane, and the plateau of adsorption isotherm for acetone in mixture is weaker than its value in the pure case. A low ratio of heptane in the mixture strongly reduces the adsorption of acetone. The influence of a little heptane on acetone adsorption is greater than the influence of a little acetone on heptane adsorption.

For an equimolar mixture (Figure 6c), the co-adsorption isotherm is under of the plateau of the adsorption isotherm of heptane alone. For this equimolar mixture, competitive co adsorption is unfavourable to both components.

These results emphasise the importance of polarity: even if acetone is more volatile than heptane, acetone is adsorbed during competitive adsorption. Electrostatic interactions between acetone and specific site (cation) are stronger than the electrostatic interaction between heptane and specific sites, so desorption of acetone from this type of site does not occur in presence of heptane.

Selectivity

The adsorption selectivity of acetone with respect to heptane 11ace/11ep is defined by the equation:

1 .

 1 used technique for high levels of contamination. The stream is oxidised at high temperature to destroy ail the pollutants completely. Its ease of application and the possibility of treating mixed gases have made this technique successful. Its main disadvantage is the production of toxic exhaust when gases contain sulphur, chlorine or fluorine. Biofiltration Here, VOCs are oxidised into CO 2 , H 2 O and sait by microorganisms fixed on an inert porous support. Biofilters are very often beds of compost and tree bark, peat or soi!. Its low cost and good efficiency at weak concentrations are the main advantages, high pressure drop and the size of installations being the main disadvantages. Recovery techniques In absorption towers, VOCs are transferred from the gas Stream into a liquid, which is often water. It is widely used in deodorization problems or for soluble VOCs, where its efficiency is good. It can also be the first stage of a hybrid process [2). The main disadvantage of this technique is the production of waste water. Condensation is used for high concentration emissions and low flow rates. The waste air is cooled to a temperature below the dew point of the mixture. The energy cost can be very prohibitive, especially for the low concentrations. High efficiency can be reached. Adsorption techniques are one of the most common methods used in solvent activities with low loaded VOCs emission in a large flow. The gaseous stream is passed through a fixed bed where VOCs are captured. The recovery of the collected vapours can be done by desorption with steam, hot air or inert gas. Activated carbon is commonly used today [3-5) but zeolites are an interesting alternative porous material for the removal of VOCs: as molecular sieves, they are not flammable, and their adsorption capacity remains good under humid conditions. Zeolites are c ry stalline aluminosilicates consisting of a breakthrough curves and adsorption isotherms are studied for three-dimensional network. The nature of a zeolite is single or mixed VOCs. The composition of the VOC mixture is characterised by its aluminium content in a lattice framework. also studied to estimate its influence on efficiency. Zeolites having a high SiO 2 / AIP 3 ratio show hydrophobie properties and can be used for the adsorption of VOCs from PILOT PLANT EXPERIMENTAL polluted air. CONDITIONS The aim of this paper is to study VOC co-adsorption on hydrophobie zeolites with a view to cleaning air. After a literature review and a presentation of the experimental conditions, the experimèntal results will be discussed. LITERATURE Ve ry few publications exist on co-adsorption of VOCs on zeolites. According to most of those found, the sieving effect is the main parameter that operates on selectivity. If the size of a molecule is greater than the aperture of the pores, adsorption does not occur [6]. For adsorbed components on a zeolite, adsorption capacities and selectivity depends on several parameters : the boiling temperature, the polarity, and the diffusion coefficient (6). The main difference between zeolites and activated carbon is the presence of cation exchange in the framework of the zeolites. These cations are specific adsorption sites that can adsorb polar components because of high electrostatie interaction between local electrostatie field (cation) and the permanent dipole moment of the co�nt {6,7). So, zeolites have two types of adsorption sites: local specific adsorption sites, represented by eations, and diffuse non specific adsorption sites, created by the framework as a whole. In the case of adsorption of chlorinated VOCs on hydrophobie faujasite, Clausse (8) concludes that, for components which present the same chemical function, the Jess volatile component is the more adsorbed and the more volatile component is the Jess adsorbed. Selectivity depends only on boiling temperature. As for co-adsorption on activated The experimental set up for adsorption is shown in Figure It consisted of an adsorption column with a fixed bed of hydrophobie synthetie zeolites. Adsorbents were in the form of pellets (HISIV from UOP) of which the main characteristies were pore diameter size dp= 6xl0• 1 0 m, specific area: 1180 m 2 g• 1 ; pellet size : 1x10-3 m diameter and 5.10-3 m height. The characteristies of the adsorption column were : 0.5 m height; 0.052 m diameter and 0.15 or 0.2 m fixed bed height. The air /VOC mixture was obtained by the use of two separate bubblers for each compound. The outlet gases were diluted with d ry air and injected at the bottom of the column. Two different VOCs were used : heptane (99%) and acetone (99.9 %). Their main physieal characteristies are given in Table 1. The desired mixture ratio was obtained by changing the flow rate of each bubbler. Three mixtures of acetone/heptane were studied, with the following molar compositions: et al. (9) studied a binary gas of m-xylene and p xylene, whieh present a difference only in their shapes. The influence of the initial composition of the gas mixture on adsorption was tested. Their results show that the selectivity of Y zeolite depends on the filling and the composition of the mixture. This is different from activated carbon adsorption. So, the aim of the present work is to show the behaviour of a commercial hydrophobie zeolite during the adsorption of solvents that are used in industry. Heptane and acetone were chosen because of their different volatilities and polarities. According to Kiselev'sclassification [10, 11] acetone is in group B and heptane in group A; for adsorbents, zeolite is of type II. That means that acetone has specific and non specific interactions with the zeolite, heptane has only non specific interactions with the zeolite. Usually, the plot of breakthrough curves is used to focus the selectivity of an adsorbent toward an adsorbate gas mixture (12, 13], and

6 Figure 1 .

 61 Figure 1. Pilot plants.1 air inlet -2,3 flowmeter -4 thermostat tank -5 bubblers -6 inlet gas -7 gas outlet sample 8 adsorption bed -9 pressure gauge -10 outlet gas.

ForFigure 2 .Figure 3 . 3 •

 233 Figure 2. Breakthrough curve for an equimolar mixture of heptane / acetone Co(ace)=0.28 moI.m-3 Co (hep )=0.26 mol m-3 . ■ C hep/Co hep ecaoe/Coaoe

2 3 4 Figure 4 . 3 •

 2443 Figure 4. Breakthrough curves for heptane (a) and acetone (b), for a same component in the same mixture for different concentrations. (a) heptane in a heptane rich mixture: ■ Co(hep) = 0.032 mol m• 3 • Co(hep) = 0.38 mol m• 3 A Co(hep) = 0.077 mol m• 3 (b) acetone in a heptane rich mixture: ■ Co (ace) = 0.038 mol m• 3 • Co (ace) = 0.028 mol m• 3 Â Co (ace) = 0.022 mol m• 3

Figure 5 .

 5 Figure 5. Partial isotherm (298 K) for heptane (a) and acetone (b) in mixture in different ratios. (a) heptane; (b) acetone 0 pure heptane 0 pure acetone • acetone-rich, heptane-poor mixture ■ acetone-poor, heptane-rich mixture  equimolar mixture

Figure 6 .Figure 7 .Figure 8 .

 678 Figure 6. Co-adsorption isotherm (298 K) and partial isotherm for heptane and acetone in three different mixtures, heptane rich mixture (a), acetone-rich mixture (b), equimolar mixture (c). 0 pure acetone D pure heptane .& global (acetone+heptane) • acetone in mixture ■ heptane in mixture


Nadhir Lebaal
email: lebaal_nadir@yahoo.fr

Fabrice Schmidt

Stephan Puissant

Design and optimization of three-dimensional extrusion dies, using constraint optimization algorithm

Keywords:

Finite element analysis

Balancing the distribution of flow through a die to achieve a uniform velocity distribution is the primary objective and one of the most difficult tasks of extrusion die design. If the manifold in a Coat-hanger die is not properly designed, the exit velocity distribution may be not uniform; this can affect the thickness across the width of the die. Yet, no procedure is known to optimize the coat hanger die with respect to an even velocity profile at the exit. While optimizing the exit velocity distribution, the constraint optimization algorithm used in this work enforced a limit on the maximum allowable pressure drop in the die; according to this constraint we can control the pressure in the die. The computational approach incorporates three-dimensional finite element simulations software Rem3D ᭨ and includes an optimiza- tion algorithm based on the global response surfaces with the Kriging interpolation and SQP algorithm within an adaptive strategy of the search space to allow the location of the global optimum with a fast convergence. The optimization results which represent the best die design are presented according to the imposed constraint on the pressure.

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. 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.

Often, the number of the 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. 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 coat hanger 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 [START_REF] Puissant | Two-dimensional multilayer coextrusion flow in a flat coat-hanger die. Part I: modeling[END_REF]. Smith et al. [START_REF] Smith | Optimal design for polymer extrusion. Part I: sensitivity analysis for nonlinear 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 shape of the extrusion die for two different materials at various temperatures, used and compared two optimization algorithms with constraint based on SQP and sequential linear programming (SLP). 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 random fluctuation function was introduced to enforce a limit on the maximum allowable pressure drop in the die.

The optimization algorithm must be carefully chosen when one single analysis using three-dimensional software requires several hours of CPU time. Non-deterministic or stochastic methods such as Monte Carlo method and genetic algorithm [START_REF] Michaeli | Approach of automatic extrusion die optimisation[END_REF] can obtain global minimum but they need a lot of evaluations for the functions to converge. Gradient methods [START_REF] Smith | Optimal design for polymer extrusion. Part I: sensitivity analysis for nonlinear systems[END_REF][START_REF] Smith | Optimal design for polymer extrusion. Part II: sensitivity analysis for weakly-coupled nonlinear steady-state systems[END_REF][START_REF] Smith | An optimisation-based approach to compute sheeting die designs for multiple operating conditions[END_REF][START_REF] Hurez | A new design procedure for profile die[END_REF][START_REF] Sun | Optimization of a flat die geometry[END_REF][START_REF] Ettinger | Parameterization and optimisation strategies for the automated design of U PVC profile extrusion dies[END_REF] require the computations of the gradients of the functions; the computation of gradients by finite difference is time consuming and depends on the perturbed parameters. For the above reasons we decided to chose a response surface method (RSM) [START_REF] Lebaal | Rheological parameters identification using in-situ experimental data of a flat die extrusion[END_REF].

The ultimate goal of this work is to optimize the coat hanger sheet die geometry (Fig. 1) in a way that a uniform velocity distribution is obtained at the die exit with an imposed nonlinear constraint so that the pressure loss in the die must decrease compared to the initial die.

For this end, we developed an automatic optimization algorithm based on SQP algorithm and a RSM together with Kriging interpolation and several strategies to permit to obtain a precise global optimum with a fast convergence. A preliminary study based on Taguchi's design of experiments method [START_REF] Chen | Optimization of the coat hanger manifold via computer simulation and orthogonal array method[END_REF] was conducted, in order to identify the most sensitive design variables. To compute a 3D flow in extrusion dies we used FEM software (REM3D ᭨ ) [START_REF] Silva | Viscoelastic compressible flow and applications in 3D injection molding simulation[END_REF]. This software takes into account strain rate and temperature dependence.

Modeling and simulation

The extrusion simulation is carried out using the 3D computation software by finite elements REM3D ᭨ [START_REF] Silva | Viscoelastic compressible flow and applications in 3D injection molding simulation[END_REF].

The flow equations are derived from the Navier-Stokes incompressible equations. A mixed finite element method for incompressible viscous flow is used. The flow solver uses tetrahedral elements with a linear continuous interpolation of both the pressure and the velocity and a bubble enrichment of velocity.

The mass, momentum and energy conservation equations, are used to follow the material behavior, from which the velocity, pressure and temperature fields are determined.

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ∇(2 ( ˙ )˙ (v)) -∇p = 0 ∇ • ⃗ v = 0 ĉ dT dt = -∇ • q + : ˙ (v) (1)
The behaviors laws used in Rem3D ᭨ give an expression of the viscosity in function of the shear rate and temperature. In this paper, the geometry of a flat die is optimized for an acrylonitrile butadiene styrene (ABS, Astalac EPC 10000). 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

= 0 (T) 1 + 0 (T) ˙ s m-1/ (2)
In this model, a, m, are material constants, whereas 0 (T) establishes the thermal dependency, given by the WLF model:

0 (T) = 0 (T ref ) exp A 1 (T ref -T s ) A 2 + (T ref -T s ) - A 1 (T -T s ) A 2 + (T -T s ) (3)
where A 1 , A 2 are material constants, and T ref is the references temperature.

A flow of 50 000 mm 3 /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.

Formulation of the optimization problem

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( ) = E E 0 such that g = P -( * P 0 ) ( * P 0 ) 0 (4)
where (J), the normalized objective function, the velocity dispersion (E), is defined as 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

E = ⎛ ⎝ 1 N N i=1 |v i -v| v ⎞ ⎠ (5)
v = 1 N N i=1 v i (6)
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 = 1; the pressure should not increase compared to the initial pressure. If < 1; the pressure must be even lower compared to the initial pressure.

Design variables

The design of a flat die is based on the values of various geometrical parameters; the optimization variables (A, B, C and D) 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 and the fourth variable D represents the thickness of the relaxation zone (Fig. 1). During the optimization process the variables A, B, C and D vary in a limited field. Each variable has its geometrical limitations which are, respectively, 30 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 (N i ), belonging to both surfaces S1, S2 and S3 (Fig. 2). For this, we use the diffuse approximation [START_REF] Ettinger | Parameterization and optimisation strategies for the automated design of U PVC profile extrusion dies[END_REF] with a linear interpolation.

A
To change the geometry of the channel repartition we define the new surface (S1 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 S1, so that they will be relocated on to the new S1 k surface which represents the new geometry. For S2 and S3 surface, a simple translation of the following nodes is applied there. The same procedure is applied on the other side of the die.

Optimization procedure

The objective and constraint function (Eq. ( 4)) are implicit compared to the optimization parameters and its evaluation need a complete numerical analysis of the extrusion process, which requires a large computation times. To need the optimum parameters with low cost and with a good accuracy, the RSM is adopted and coupled with an auto adaptive strategy of the research space Fig. 3. 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 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 (Fig. 3). 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 = 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 =10 -3 .

Kriging interpolation

The Kriging interpolation [START_REF] Kaymaz | Application of Kriging method to structural reliability problems[END_REF][START_REF] Li | Development of a novel meshless local Kriging (LoKriging) method for structural dynamic analysis[END_REF], makes it possible to present the complexes function effectively. This method is applied in our work to represent the response surface in an explicit form, according to the variables of optimization. The approximate explicit relationship of the objective and constraint function, can be expressed as follows:

J(x) = pT (x) â + Z(x) ( 7)
with p(x) = [ p1 (x), . . . , pk (x)] T , where k denotes the number of the basis function in regression model, â = [ â1 , ... , âk ] T 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 pT (x) â in Eq. ( 7) 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. ( 7) is a correction of the global model. It is used to model the deviation from pT (x)a so that the whole model interpolates response data from the function The construction of a Kriging model can be explained as follows:

The output responses from the function are given as

F(x) = {f 1 (x), f 2 (x), . . . f n (x)} ( 8
)
From these outputs the unknown parameters â can be estimated:

â = ( PT R -1 P) -1 PT R -1 F ( 9
)
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 = ⎡ ⎢ ⎣ R(x 1 , x 1 ) • • • R(x 1 , x n ) . . . . . . . . . R(x n , x 1 ) • • • R(x n , x n ) ⎤ ⎥ ⎦ + ⎡ ⎢ ⎢ ⎢ ⎣ w(x -x 1 ) 0 • • • 0 0 w(x -x 2 ) • • • 0 . . . . . . . . . . . . 0 0 • • • w(x -x n ) ⎤ ⎥ ⎥ ⎥ ⎦ ( 10
)
R ij = (|x i -x j |) (11)
A weight function of Gaussian type with a circular support is adopted for the Kriging interpolation. It takes the form

w i (x) = ⎧ ⎪ ⎨ ⎪ ⎩ 1 - e -(d i /c) 2 -e -(rw/c) 2 1 -e -(rw/c) 2 si d i r w 1 si d i r w ( 12
)
where 2 is the distance from a discrete node x i to a sampling point x in the domain of support with radius r w , and c is the dilation parameter. c = r w /4 is used in computation.

d i = n J=1 (x J -x J (i) )
The second part in Eq. ( 7) is in fact an interpolation of the residuals of the regression model pT (x) â. Thus, all response data will be exactly predicted; is given as

Z(x) = r T (x) ( 13
)
where r T (x) = {R(x, x 1 ), . . . , R(x, x n )}. The parameter is defined as follows:

= R -1 (F -P â) (14)

Result and discussion

Design sensitivity analysis

A Taguchi's design of experiments method [START_REF] Chen | Optimization of the coat hanger manifold via computer simulation and orthogonal array method[END_REF] is used to investigate the sensitive analysis of the design variables. This preliminary study was conducted, by considering four geometrical parameters (A, B, C and D), in order to identify the most sensitive ones on the objective and the constraint function.

For this, we have calculated the effect of each variable over all the field of research by studying the contribution of each variable on the variation of the response "velocities distribution and pressure". 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.

We note in Fig. 4 that the variable with less influence on the objective function is the variable D which represents the thickness of the relaxation zone. On the other hand, the variable C has the greatest effect on the velocity distribution at the die exit.

We observe in the seam figure that if the variables "A, B, and C" are in the lowest level (1), the global relative variation is more important which implies a bad velocity distribution at the die exit. On the other hand, if the variables are in the highest level (2), the global relative variation of velocity at the die exit is less important compared to their lowest level.

For the variable "C", it is even weaker at the average level, and for the variable "D", it is noted that it has a very weak effect over the velocity distribution at the die exit.

From Fig. 5 we note that the effect of the variables B and D on the constraint function (pressure) is weak. We note that the more the variables "A, C and D" are in the lower level the more the pressure increases. On the other hand, the pressure is lower at their highest level. For the variable B, the variation of the pressure is very weak. Since the prime objective in the design of the extrusion dies is to obtain at the exit a perfectly homogeneous product, it was noted that the variable D has the least important effect on the velocity distribution at the die exit. This study enables us to limit the number of variables and three sensitive geometrical parameters are selected and considered as design variables to be optimized (A, B and C).

Optimization results

The optimization strategy described previously is applied to optimize the die geometry in order to obtain a homogenous exit velocities distribution. According to the imposed constraint on the pressure, two optimization cases are formulated:

Case 1: A constraint on the pressure is imposed for which the pressure in the optimal die must be lower or equals 0.5 * P 0 ( = 0.5).

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

A representative convergence history during an optimization process on the objective functions (J) is represented in Fig. 6. After three iterations, the objective function fell below the cutoff of 10 -3 and the simulation stopped.

Given the fixed geometry constraints of the die, the rheological parameters of the ABS resin, and the process conditions, the optimal solution is obtained starting from the second iteration for both cases. The objective function is then reduced to 95% of its initial value for the first case ( = 0.5) and 91% for the second case ( = 0.35). The optimization strategy clearly showed its capacity to obtain an optimal solution with a very fast convergence. A summary of the optimization results obtained for the two cases, is reported in Table 2.

We notice that the objective function for both cases 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 case 1 offer the best minimum, because of the less severe constraint on the pressure compared to case 2.

We have represented the relative velocity dispersion obtained in the initial geometry and after optimization, for both results in Fig. 7. We see that the objective function obtained by imposing a severe constraint on the pressure ( = 0.35) is slightly higher of the objective function obtained with a less severe constraint ( = 0.5). For a qualitative comparison it is obvious that velocity distribution in the initial geometry (Fig. 8A) is not homogeneous. We notice that in the middle of the die, velocities are lower about 35% of variation compared to the average. On the other hand, in both parts close to the end of the die, we observe an increase of velocities until 50% (Fig. 7), only in the small regions near the end of this die, the exit velocity reduces to zero to satisfy the no slip condition at the walls. After optimization, we notice that velocities became homogeneous all over the width at the die exit. For both optimization results, seeing, respectively, Fig. 8B andC. The constraint on the pressure is represented on Fig. 8, in a form of (P/P 0 ) for both cases. We observe that the constraint function is weaker with a more severe constraint ( = 0.35).

This optimization problem was formulated in a way to find the best velocity distribution at the die exit while avoiding the increase of pressure compared to an imposed pressure, which is respected by both results (Fig. 9). The pressure drop in the design of the initial die is 4.42 MPa (Fig. 10A). It should be noted that the pressure decreases during the iterations of optimization for the different results. In the first case, the pressure decreases by 52% which gives a pressure of 2.14 MPa (Fig. 10B). In the results of the second case, we note that the pressure has decreased by 62% because of the imposed constraint (Fig. 9), which gives a pressure of 1.69 MPa (Fig. 10C). We notice that the gradient pressure in the initial die is not parallel, which results on a non-uniform flow. Nevertheless, in the optimized dies the gradient pressure is parallel; furthermore, it denotes a uniform flow distribution for both cases.

It is noted that more the polymer runs towards the border because of the heating by shearing, the temperature is lower in the middle of the die and higher close to the border (Fig. 11A). By comparison, the distribution of the temperatures is uniform in the optimal die (Fig. 11B andC), because of the homogeneous velocities distribution (due to the heating by shearing) at the die exit for the optimal geometries.

Conclusion

The design of polymer extrusion dies is complicated by the nonlinear relationship between the resin viscosity, shear rate and temperature. The first study based on the sensitive analysis enabled us to decrease the number of variables by analyzing their effects on the pressure and the velocity distribution at the die exit. The numerical optimization algorithm presented in this work shows its robustness as a tool for extrusion die design. The applied optimization method with the various strategies proves its ability of predicting optimal die geometry with satisfactory computing time for three-dimensional finite element calculation.

The obtained results with the two optimization examples show that our optimization strategy is well adapted to the nonlinear constraints. The use of the different constraints on the pressure allowed us to compare the various results obtained after optimization. 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 algorithm of optimization clearly showed its robustness and its capacity to obtain an optimal solution with a fast convergence.

Future research will involve experimental verification of the algorithm and the optimization of complex helicoidally dies.

3 Fig. 1 .

 31 Fig. 1. Geometry of the flat die.

Fig. 2 .

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

Fig. 3 .

 3 Fig. 3. The flowchart of the adopted optimization strategy.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Effect of the variables on the objective function.

Fig. 6 .

 6 Fig. 6. Convergence history during an optimization process on the objective functions.

Fig. 7 .

 7 Fig. 7. Velocity dispersion with initial and optimal die.

Fig. 8 .

 8 Fig. 8. Velocity distribution in initial and optimal dies.

Fig. 9 .

 9 Fig. 9. Convergence history of the constraint functions during an optimization run.

Fig. 10 .

 10 Fig. 10. Pressure distribution in initial and optimal dies.

Fig. 11 .

 11 Fig. 11. Temperature distribution in mid plane in initial and optimal dies.

Table 1

 1 Rheological parameters

	20.4	101.6	397.7	524.7

0 (Pa s) m 0 s (Pa) 716.997 0.15862 1 224063 A1 A2 (K) Ts (K) Tref (K)

Table 2

 2 Summary of the results

	Results of optimization	Initial die (A)	Result 1	Result 2
			= 0.5 (B)	= 0.35 (C)
	"J"	1	0.049	0.09
	Improvement (%) "g"	-1	95 0.48	91 0.38
	Profit of pressure (%)	-	52	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


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

Yannick Le Maoult

S Monteix

Modelling of infrared heating of thermoplastic sheet used in thermoforming process

Keywords: Thermoforming process, Polystyrene, Halogen lamps, Surface temperature measurement, Infrared camera, Control volume method

come

Modelling of infrared heating of thermoplastic sheet used in thermoforming process

Introduction

In Fig. 1, the plug-assisted thermoforming process is sketched. The heating zone is constituted of upper and lower rows of infrared emitters, which are used to warm the plastic sheet prior to forming [START_REF] Illig | Pratique du thermoformage[END_REF][START_REF] Throne | Two-Day Seminar on Thermoforming Process and Design[END_REF].

In most of the thermoforming machines, the heating step is performed using an infrared oven constituted of long waves infrared emitters (3-5 m spectral bandwidth). The goal of this study is to determine the efficiency of short waves infrared emitters (halogen lamps) for the heating step (0.76-2 m spectral bandwidth) in the thermoforming process. For the purpose, a representative polymer of the thermoforming process has been chosen, i.e. the polystyrene (PS). According to Jammet [START_REF] Jammet | Thermoformage, Techniques de l'ingénieur[END_REF], this polymer represents 49% of the entire thermoformed polymers. Then, comparisons will be made upon energetic criterion between three different emitters: long waves infrared emitters (ELSTEIN), average waves infrared emitters (TQS) and halogen lamps (PHILIPS).

Experimental set-up

Fig. 2 shows the experimental set-up, which has been developed in order to measure the surface temperature evolution when a white PS sheet is heating using an infrared lamp (one or more if necessary). The polystyrene sheet is mounted in an isolating frame. The surface dimension of the sheet is 135 mm × 100 mm and the thickness is 1.63 mm. An 880 LW AGEMA I.R. camera (8-12 m bandwidth) is used to evaluate the surface distribution of the transmitted heat flux by measuring the temperature distribution on the back surface of the sheet. The frequency of analysis is 25 frames/s and the device is connected to a 12 bit data acquisition board drive with a real time software. The experimental set-up allows for monitoring the transient surface temperature distribution. An aluminium polish sheet with an aperture (diaphragm) can be placed between the polystyrene sheet and the lamp in order to insure a constant view factor between the different lamps and the polymer sheet.

Heat source characterisation

The halogen lamps used in the infrared oven (PHILIPS 900 W) are composed of a coiled tungsten filament (length: 300 mm), contained in a quartz tubular enclosure (diameter: 10 mm) filled with a neutral gas (Argon) and coated on its back with a ceramic reflector in order to increase the heat flux received by the product. This lamp as well as the EL-STEIN and the TQS lamps are presented in Fig. 3. Spectral properties of these elements have already been measured in previous works [START_REF] Monteix | Experimental Study and Numerical Simulation of Sheet and Tubular Preform Infrared Heating, QIRT[END_REF][START_REF] Monteix | Recent Issues in Preform Radiative Heating Modelling, PPS 2001[END_REF]. Using these data, the tungsten filament and quartz tube temperatures are calculated using a net radiation method [START_REF] Petterson | Heat transfer and energy efficiency in infrared papers dryers[END_REF]:

P -ε Fil (T Fil )S Fil σT 4 Fil -2πL Q k Argon (T * ) T Fil -T Q ln(d Q /d Fil ) = 0, α Q (T Fil )ε Fil (T Fil )S Fil σT 4 Fil + 2πL Q k Argon (T * ) T Fil -T Q ln(d Q /d Fil ) -ε Q S Q σT 4 Q -hS Q (T Q -T ∞ ) = 0 ( 1
)
where P is the electrical power of the lamps; the indices Fil and Q refer, respectively, to the filament and the quartz tube, T the temperature, L the length of the lamp, d the filament or the quartz tube diameter, ε, α the emissivity and absorptivity. T * is assumed to be equal to (T Fil + T Q )/2. We have simplified this model by computing directly total parameters instead of spectral ones (as emissivity and absorptivity) with polynomial functions depending on temperatures which is easier and saves time. The temperatures of the filament and the quartz tube are computed at steady state using a Newton-Raphson method. The results of this computation are shown in Fig. 4. Additional curves are plotted on the same graph, they are related to experimental measurements: colour temperature of the filament and temperature computed with a resistivity method. At the nominal power of 1 kW, the values are given in Table 1.

As shown in Fig. 4, the discrepancy between the extreme curves is high: greater than 100 K around 2350 K. This fact is due to the lack of knowledge on the properties of this particular kind of tungsten (coiled filament). On the other hand, "colour" temperature remains difficult to measure at such level of power. In order to conclude, we have chosen the computed temperature as a "mean" temperature, lower and upper limits of the temperature, at 1 kW, will be considered as the global uncertainty.

Using experimental set-up described in [START_REF] Monteix | Experimental Study and Numerical Simulation of Sheet and Tubular Preform Infrared Heating, QIRT[END_REF][START_REF] Monteix | Recent Issues in Preform Radiative Heating Modelling, PPS 2001[END_REF], the efficiency factor of the ceramic reflector has been measured:

k r (θ, T r ) = Φ r (θ, T r ) Φ w (θ, T r ) (2)
The index r and w refer, respectively, to the lamp with a ceramic reflector and without reflector. The angle θ covered by a detector around the heater and T r the temperature of the reflector. This value is obtained using a thermopile detector with a large spectral bandwidth [0.25, 26] m. The output signal delivered by the thermopile is related to the flux emitted by the lamp, so the ratio of the two fluxes measured by the detector and introduced in (2) gives the efficiency of the coating. We have checked that this efficiency is unaffected by the variation of the temperature of the filament and that the emission of the lamp is quasi isotropic. Therefore, we have obtained

k r = 1.36 ± 0.02 (3)

PS spectral properties characterisation

Spectra data will be used as a database in the numerical simulation (see Section 3) as well as for the surface temperature distribution measurements. In the literature, we found few papers referring to polystyrene spectral properties measurements [START_REF] Esser | Infrared radiation in the processing of plastics: precise adjustment-the key to productivity[END_REF][START_REF] Miyanaga | Analysis of Infrared Radiation Heating of Plastics[END_REF]. 1.63 mm white PS sheets were processed by extrusion moulding. In addition, the previous thick sheets were polished in order to obtain samples of 0.07 mm thickness. Then, we have proceeded to spectral properties measurements (transmitivity, reflectivity) using a Perkin-Elmer FT-IR spectrometer. We have found that the polymer is opaque for wavelengths in the range of 8-12 m with an integrated emissivity of 0.95. This bandwidth corresponds to the one of the infrared camera. In Fig. 5, the PS absorbtivity α λ is plotted versus wavelength both for the 1.63 mm sheet (a) and the 0.07 mm thin film (b).

Numerical modelling

3D control-volume software, called PLASTIRAD, has been developed in order to compute heat transfer during the infrared-heating step. The sheet is meshed using cubic or hexahedral elements so-called control volumes [START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF]. The temperature balance equation including radiant heat transfer (thermal power absorbed by the semitransparent PS sheet

= - t Γ e (⃗ q c • ⃗ n) dΓ dt + t Γ e (⃗ q r • ⃗ n) dΓ dt (4)
where Ω e is the control volume, q c and q r are the Fourier's conductive and the radiated fluxes, respectively. The unknown temperatures are computed at the cell centres of the elements. The specific lamp geometry is also considered, taken into account from a view factor computation.

In our model, we assume a constant filament temperature and then a uniform source temperature. Thus the amount of radiation, reaching each front elements of the irradiated sheet, is performed with the contour method using Stokes theorem [START_REF] Rammohan | Efficient evaluation of diffuse view factors for radiation[END_REF]:

F H→p = 1 2πS H Γ H Γ rp ln r • d⃗ r H • d⃗ r r p ( 5
)
where the subscripts H refers to the heater and p to the sheet front element under consideration, r is the distance between two elemental linear elements, issued from meshing the entire heater boundary Γ H and Γ r p for each surface of control element's boundary reaching the irradiated sheet's face. This contour method integration was chosen between different numerical view factor methods (Monte Carlo or geometrical method). This method allows a good accuracy with a very low computational time. An efficient numerical method based on Gaussian quadrature has been applied in order to perform the contour integration. The computation of the amount of incident radiation, for each surface element, is given in Eq. ( 6):

Φ(r p ) = k r F H→r p S H S r p λ ε λ (T Fil )πL 0 λ (T Fil ) dλ (6)
where ε λ is the spectral tungsten emissivity, T Fil the filament temperature, L 0 λ the blackbody intensity, λ a given wavelength between [0.2, 10] m and k r the efficiency factor of the ceramic reflector. In a first approximation, internal radiant heat transfer is assumed mono-dimensional across the sheet thickness. The PS bulk temperature (≤450 K) is very low in comparison to the source temperature (T Fil > 2000 K). Therefore, the assumption of cold material is convenient. This leads to express the transmitted flux across the thickness of the material, for an homogenous and amorphous PS, using the Beer-Lambert's law, if the coordinate y describe the thickness, we have, in a simple way, for the flux

Φ λ (y) = Φ λ (y = 0) e -k λ y (7)
where k r is the spectral absorption coefficient of our material.

Discussion

Model assessment

First of all, we proceed to a heating step of a white PS sheet during 70 s using one halogen lamp. The distance between the lamp and the sheet is 50 mm. It is to be noted that in this case, the diaphragm is not taken into account. The dimensions of the sheet are the same as in Section 2. The sheet is meshed using 2250 linear cubic elements (15 × 15 × 10). Fig. 6 shows characteristic back surface temperature distribution obtained by numerical simulation (a) and by experimental measurement (b). We obtained a qualitative agreement. More interesting is the temperature variation along the width and the height of the sheet from a cross-point located on the centre of the sheet (see Fig. 6a). The maximum relative error for the temperature distribution along width (Fig. 7a) is 2.5%. In this case, the agreement is fair. Moreover, the relative error in the case of the temperature distribution along height (Fig. 7b) increases from 4% at the sheet centre up to 35% at the sheet border. This discrepancy is due to the mono-dimensional absorption assumption in our model (see Section 3).

Comparison between different emitters

In this section, we compare three different emitters: a long wave infrared emitter (ELSTEIN 650 W), an average wave infrared emitter (TQS 600 W) and a halogen lamps (PHILIPS 900 W) previously introduced in Section 2.

We proceed to a heating step of 290 s for each emitter using the experimental set-up described in Section 2. This is the duration for the less efficient infrared emitter to achieve a maximum temperature of 74 • C. The aluminium polish sheet with an aperture (diaphragm) is placed between the polystyrene sheet and the lamp in order to avoid geometric effect. Then, the temperature evolution is recorded versus time at the sheet centre using the infrared camera (see Fig. 8). Then, we calculate a simple energetic criterion

E L = P L × t L ( 8
)
where E L is the approximated energy consumed by each lamp, t L the time required for each lamp to achieve the less efficient infrared emitter maximum temperature, and P L the electrical power of each lamp. If we applied Eq. ( 8) to the curves of Fig. 8a, we obtain the values refereed in Table 2.

According to the criterion (Eq. ( 8)), the energy efficiency of the halogen lamp is better than the long waves infrared emitters (i.e. ELSTEIN lamp). In addition, the PHILIPS lamp efficiency should be improved in comparison to that of the TQS lamp. For that, an aluminium reflector (IRZ 1000) has been placed in the back of the halogen lamp in order to increase radiant heat flux. In Fig. 8b, the temperature evolutions of the different lamps are plotted versus time. Again, we calculate the energy consumed by each lamp (Table 3).

Due to the presence of the aluminium reflector, the halogen lamp has now the best energy efficiency (i.e. lower energy consumed).

Conclusion

An experimental set-up has been developed in order to measure the back surface temperature distribution of a polystyrene sheet using an infrared camera and has permitted to compare the energy efficiency of three different emitters (TQS, ELSTEIN, PHILIPS). The comparisons, based upon a simple energetic criterion, allow demonstrating that the halogen lamp has the best energy efficiency. These preliminary results will be further improved and extended to other polymers, and particularly to polypropylene (PP).

A 3D numerical model using control volume method has been developed to simulate the heating stage and then validated by comparisons with experimental data. This numerical model has to be improved in order to take into account the transmitted flux across the thickness without geometric assumptions (i.e. Beer-Lambert mono-dimensional law).

Fig. 1 .

 1 Fig. 1. Plug-assisted thermoforming process.

Fig. 2 .

 2 Fig. 2. Experimental set-up overview. Side view of the infrared oven.

Fig. 3 .

 3 Fig. 3. PHILIPS (with and without reflector), ELSTEIN and TQS lamps.

Fig. 4 .

 4 Fig. 4. Filament temperature vs. electrical power.

Fig. 5 .

 5 Fig. 5. Absorptivity vs. wavelength.

Fig. 6 .

 6 Fig. 6. Back surface temperature distribution.

Fig. 7 .

 7 Fig. 7. Numerical and experimental temperature distribution along width (a) and along height (b).

Fig. 8 .

 8 Fig. 8. Recorded temperature vs. time at the sheet centre.

Table 1

 1

	Filament and quartz temperatures
	Filament temperature (K)	Quartz temperature (K)
	2360	730

Table 2

 2

	Energy consumed by each lamp
		TQS	PHILIPS	ELSTEIN
	E L (kW)	54	144	188.5
	Table 3
	Energy consumed by each lamp
		TQS	PHILIPS	ELSTEIN
	E L (kW)	54	49.5	188.5

Acknowledgements

PHILIPS LIGHTING and EDF-RENARDIERE Companies support this work. Special thanks go to R. Anbari for his help on the experiments.