V Velay

L Robert

F Schmidt

S Hmida

T Vallet

Behaviour model identification based on inverse modeling and using Optical Full Field Measurements (OFFM) : application on rubber and steel

Keywords: Full field measurements, Forward problem, Behaviour identification PACS: 62.20.-x

Biaxial properties of materials (polymer or steel) used in many industrial processes are often difficult to measure. However, these properties are useful for the numerical simulations of plastic-processing operations like blow moulding or thermoforming for polymers and superplastic forming or single point incremental forming for steels. Today, Optical Full Field Measurements (OFFM) are promising tools for experimental analysis of materials. Indeed, they are able to provide a very large amount of data (displacement or strain) spatially distributed. In this paper, a mixed numerical and experimental investigation is proposed in order to identify multi-axial constitutive behaviour models. The procedure is applied on two different materials commonly used in forming processes: polymer (rubber in this first approach) and steel. Experimental tests are performed on various rubber and steel structural specimens (notched and open-hole plate samples) in order to generate heterogeneous displacement field. Two different behaviour models are considered, on the one hand, a Money-Rivlin hyperelastic law is investigated to describe high levels of strain induced in tensile test performed on rubber openhole specimen. On the other hand, Ramberg-Osgood law allows to reproduce elasto-plastic behaviour of steel on notched specimen. Each model identification is based on a same Finite Element Model Updated (FEMU) procedure which consists in comparing results provided by numerical simulation (ABAQUS T M ) with full field measurements obtained by DISC (Digital Image Stereo-Correlation) technique (Vic 3D ). In both cases, a cost function is defined minimizing in a least square sense the difference between experimental and numerical displacement fields induced into the structural specimen. In this paper, experimental tensile tests performed on rubber and steel specimen will be described, displacement and strain field measurements obtained by DISC technique will be analyzed. Then, the different steps of the methodology will be detailed. Displacement field calculated by numerical simulation will be available at each node of the FE meshing and will be properly associated to the experimental values measured by the optical full field method. Afterwards, these information will be useful to define the cost functions. Moreover, boundary conditions applied in the simulation, investigation area selected to perform identification, rigid body motion included in the experimental data and cost function forms will be discussed. Finally, results provided by the procedure will be presented for both investigated materials.

INTRODUCTION

The aim of this investigation is to present an original methodology able to determine multiaxial material properties. The approach can be applied to identify behaviour models of steels or polymers commonly used in many industrial forming processes like blow moulding or single point incremental forming. The paper is based on a mixed experimental and numerical approach using Optical Full Field Measurements (OFFM) and Finite Element model results. Today, this method is more and more applied in experimental mechanics, it allows to measure kinematic fields (displacement or strain) and provides extended possibilities to identify multiaxial behaviour models. Several recent investigations using different experimental full field methods and various identification procedures have been developed [START_REF] Meuwissen | An inverse method for the mechanical characterization of metals[END_REF][START_REF] Bonnet | [END_REF]. In this investigation a Finite Element Model Updated (FEMU) method is considered [START_REF] Haddadi | Identification of the parameters of Swift law using strain field measurements[END_REF][START_REF] Giton | [END_REF][START_REF] Lecompte | Elastic orthotropic parameter identification by inverse modeling of biaxial tests using Digital Image Correlation[END_REF]. It consists in comparing experi-mental measurements obtained by DISC (Digital Stereo-Correlation) technique (VIC 3D ) with numerical simulation (ABAQUS T M ). Two different kind of tensile tests are investigated on two structural specimens and two materials (rubber and steel). The aim is to induce a multiaxial loading within the specimen, thus an open hole specimen is defined to investigate rubber behaviour whereas a notched specimen is investigated for steel behaviour.

EXPERIMENTAL INVESTIGATION

Both of tensile tests performed on rubber and steel have been carried out with a servo-electric testing machine and controller connected to a computer. It allows to measure the displacement rate and the global strength induced into the specimen. Two CDD cameras are located in front of the testing machine in order to measure dis-placement field during the tests by DISC technique (Figure 1).

FIGURE 1. DISC set up and tensile testing machine

In each case, the investigation area only considers a part of the structural specimen and can not include the whole sample. This experimental method allows to provide displacement or strain field components induced within the investigation area by the global loading. Each displacement or strain map corresponds to a time increment of the tensile test. Figure 2 presents global tensile curve obtained for the steel specimen. In the following part, three images will be selected for the identification process corresponding to three time increments (100s, 200s and 300s). Similar assumption will be done for the rubber specimen. For instance, figure 3 shows axial and shear strain components within the steel notched specimen at a time of 300s. Kinematic field induces a rigid body motion will have to be taken into account into the identification process.

IDENTIFICATION METHODOLOGY Strategy

Identification strategy can be divided into several stages. First, some preliminary works have to be performed on the experimental files in order to keep only the correlated solutions and to avoid mismatched informations. Then, locations where experimental data are measured have to be projected onto the finite element mesh. Figure 4 presents for the rubber open-hole specimen a comparison between nodes considered in experimental and numerical meshes. Only an investigation area is modelized in the numerical simulation. Experimental and calculated displacements will be considered to evaluate the cost function. Moreover, experimental displacements measured at the upper and lower parts of the specimen will be introduced into the simulation to define the boundary con-ditions. Thus, rigid body motions will be taken into account in the simulation. Finite element simulation will be performed with ABAQUS T M , a script allows to extract node displacements and global strength necessary to the analysis. The main program automatically runs the finite element calculation, extracts calculated and experimental data for the cost function evaluation and updates behaviour model parameters to be identified. Cost function minimization is performed with modified Nelder-Mead Simplex method implemented in Matlab .

The cost function depends on experimental and calculated displacements and strength. Its formulation includes a spatial (nodes) and time informations (images).

J(P) = β . ||u exp -u sim || ||u exp || + (1 -β ). ||F exp -F sim || ||F exp ||
where:

||u exp -u sim || = ∑ N n ,N i k,i=1 u sim 1 -u exp 1 2 k,i + u sim 2 -u exp 2 2 k,i
and:

||F exp -F sim || = ∑ N i i=1 F sim -F exp 2 i
N n is the number of nodes considered for the analysis, N i the number of images obtained at various time increments, u i , i = 1, 2 the displacement components, F the global strength in the loading direction, P the parameter vector to be identified and β the weight parameter.

Figure 5a shows the calculated displacement magnitude field corresponding to the first image (t=100s). Applying the experimental displacement field as a boundary conditions (upper and lower parts of the specimen), rigid body displacement fields are naturally taken into account into the calculation. Figure 5b illustrates reaction forces induced by the previous boundary conditions applied. They are considered for cost function evaluation.

Behaviour modeling

Two relevant laws are carried out in order to investigate rubber and steel behaviour. Indeed, a Ramberg-Osgood law is considered to describe steel behaviour whereas a Mooney-Rivlin hyperelastic law is used to reproduce behaviour of rubber.

steel

A Ramberg-Osgood model is selected to reproduce the steel behaviour. This law is available into ABAQUS T M software. Basic one-dimensional model is:

Eε = E(ε e + ε p ) = σ + E K |σ | K 1 M -1
σ E is the Young Modulus, ε the total strain, ε e and ε p the elastic and plastic parts, σ the stress induced, K and M the parameters of the non linear term defining the hardening. Multiaxial formulation includes Hooke linear elastic relation used to generalize the first term of the previous equation and a nonlinear term to define plasticity through the use of the Mises stress potential and associated flow law.

Eε = (1 + ν)σ ′ -(1 -2ν)pI + 3 2 E K σ eq K 1 M σ ′
Elastic properties E and ν are assumed to be known and only both parameters of the nonlinear term K and M will be identified.

In order to assess the material parameters of the forward problem to be identified, various direct simulations have been performed in both cases (steel and rubber). Figure 6 presents evolution of the cost function versus the two parameters K and M of the Ramberg-Osgood law.

A global minimum seems to be reached around a value of 5 for M and 500 for K.

Rubber

Mooney-Rivlin hyperelastic model is known to provide a good description of rubber behaviour. Moreover, it only requires two material parameters C 01 and C 10 to be identified. Energy function can be formulated as follow:

Results and discussion

For both cases (steel and rubber), preliminary tensile tests performed on flat specimens allows to assess a first set of material parameters and to initialize the identification procedure.

Steel

Several minimizations are performed in order to investigate the influence of the weight parameter β . A value of 0.9 confers a more important weight on the displacement field. On the contrary, a value of 0.2 allows to maximize the effect of the force induced. Figure 8 gives evolution of the Ramberg-Osgood parameters versus the number of iterations and figure 9 the evolution of the cost function for a weight parameter equals to 0.9. A good correlation on both displacement fields and reaction forces is obtained. Figure 10 presents axial (a) and shear (b) strain components for a time of 300s. Results can be compared with those provided by DISC technique (Figure 3). They are in a good agreement with experiment.

Table 1 gives presents final results.

Rubber

In this case, DISC technique is more difficult to carry out due to the very large displacement induced in tensile tests. Presently, experimental results allow to assess the heterogeneous deformation within the specimen. But, these measurements can not yet be used to determine behaviour parameters by an inverse approach.

CONCLUSIONS

In this investigation, an efficient methodology was carried out in order to identify multiaxial behaviour laws within various test conditions and different materials. It is based on a forward problem using optical full field measurements and finite element simulations. For that purpose, a Finite Element Model Updated method was developed. A cost function was formulated in term of displacement fields and reaction forces and minimized by a modified Nelder-Mead Simplex algorithm. Two different applications were considered. On the one hand, identification of elastoplastic behaviour law of steel was successfully performed. Strain component fields provided by numerical simulation are in a good agreement with those obtained from an experimental point of view. On the other hand, the inverse approach was applied on a hyperelastic behaviour of rubber. In this case, some experimental problems occur. Test conditions have to be improved in order to be able to perform a complete multiaxial identification.

FIGURE 2 .

 2 FIGURE 2. Tensile test performed on notched steel specimen

FIGURE 3 .

 3 FIGURE 3. a-Axial strain field; b-Shear strain field measured by DISC Technique

FIGURE 4 .

 4 FIGURE 4. Finite element mesh projection onto experimental data

FIGURE 5 .

 5 FIGURE 5. a-Displacement magnitude field induced in the steel notched specimen, b-Reaction forces induced by the boundary conditions

FIGURE 6 .

 6 FIGURE 6. Evolution of the cost function versus the material parameters (Ramberg-Osgood law)

FIGURE 7 .

 7 FIGURE 7. Displacement fields provided by a multi-stage numerical simulation performed on rubber open-hole specimen)

FIGURE 8 .FIGURE 9 .

 89 FIGURE 8. Evolution of the Ramberg-Osgood parameters K and M over the number of iterations for a weight parameter of 0.9

FIGURE 10 .

 10 FIGURE 10. a-Axial strain field; b-Shear strain field calculated by Finite Element Simulation

Figure 11

 11 Figure 11 illustrates a comparison between measured and simulated shear strain field obtained with coefficients identified in uniaxial test conditions (C 10 = 0.136 and C 01 = 0.097). Strain concentration can be predicted by optical measurements represented on only one half of specimen.

FIGURE 11 .

 11 FIGURE 11. Comparison between experimental and calculated shear strain fields performed on open-hole rubber specimen

TABLE 1 .

 1 Final results provided by identification process (β = 0.9)

	E MPa (fixed) ν (fixed) K MPa	M
	166000	0.34	481	4.92


N Pirc

F Schmidt

M Mongeau

F Bugarin

F Chinesta

OPTIMIZATION OF BEM-BASED COOLING CHANNELS INJECTION MOULDING USING MODEL REDUCTION

ABSTRACT: Today, around 30% of manufactured plastic goods rely on injection moulding. The cooling time can represent more than 70% of the injection cycle. In this process, heat transfer during the cooling step has a great influence both on the quality of the final parts that are produced, and on the moulding cycle time. Models based on a full 3D finite element method renders unpractical the use of optimization of the design and placement of the cooling channel in injection moulds. We have extended the use of boundary element method (BEM) to this process. We introduce in this paper a practical methodology to optimize both the position and the shape of the cooling channels in injection moulding processes. We couple the direct computation with an optimization algorithm such as SQP (Sequential Quadratic Programming). First, we propose an implementation of the model reduction in the BEM solver. This technique permits to reduce considerably the computing time during the linear system resolution (unsteady case). Secondly, we couple it with an optimization algorithm to evaluate its potentiality. For example, we can minimize the maximal temperature on the cavity surface subject to a temperature uniformity constraint. Thirdly, we present encouraging computational results on plastic parts that show that our optimization methodology is viable. KEYWORDS: BEM, optimization, reduction model, injection moulding, SQP.

INTRODUCTION

Numerical simulations for designing injection moulds have become an important developement in injection moulding processes. The location of the cooling channels is a major element in the design of the mould because the cooling time can represents up to 70 % of the injection cycle. We need efficient numerical simulations in order to optimize the process parameters, but models based on full 3D finite element method renders unpractical the use of optimization for this design and placement of cooling channels in injection moulds. In this context, the Dual Reciprocity Method (DRM), introduced by Brebbia [START_REF] Chen | Dual reciprocity method using compactly supported radial basis functions Communications in Numerical Methods in Engineering[END_REF], is acknowledged to be one of the most effective BEM techniques for transforming domain integrals into boundary integrals. More-over, during the optimization, this method permits to compute exact gradients, thereby avoiding the N direct computations per optimization iteration that are needed by finite-difference gradient approximation (where N is the number of optimization variables). However, thermal models involved in the numerical modelling in injection moulding processes a certain number of numerical difficulties such as size meshing, long simulations, or the necessity to define a homogenized thermal conductivity. In the first part of this paper, we present the use of boundary element method (BEM) and DRM applied to unsteady heat transfer of injection moulds. The BEM software, developed at the CROMeP laboratory [START_REF] Mathey | Automatic optimization of the cooling of injection mold based on the boundary element method Materials Processing and Design: Modeling, Simulation and Applications[END_REF], was combined with an adaptive reduced modelling [START_REF] Chinesta | Alleviating mesh constraints: model reduction, parallel time integration and high resolution homogenization[END_REF]. This procedure permits to reduce considerably the computing time during the linear system resolution in unsteady prob-lem. Then, we present a practical methodology to optimize both the position and the shape of the cooling channels in injection moulding processes. We couple the direct computation with an optimization algorithm such as SQP (Sequential Quadratic Programming) [START_REF] Nocedal | Numerical optimization series in operation research[END_REF]. For the sake of simplicity, we will consider a potential problem defined in a 2D unbounded domain. The capabilities of both the reduced order modeling and the boundary element method will be outlined.

BEM AND DRM APPLIED TO POISSON EQUATION

Using BEM, only the boundary of the domain has to be meshed and internal points are explicitly excluded from the solution procedure. An interesting side effect is the considerable reduction in size of the linear system to be solved [START_REF] Pirc | BEM-based cooling optimization for 3D injection molding[END_REF]. The transient heat conduction in a homogeneous isotropic body Ω is described by the diffusion equation [START_REF] Kenig | Cooling molded parts, a rigorous analysis[END_REF], where a is the material diffusion:

∀M ∈ Ω, -→ ∇ 2 T (M, t) = 1 a ∂(M, t) ∂t (1)
We define the initial conditions and the boundary conditions as:

   T (M, t) = T 0 ∀M ∈ Γ P φ(M, t) = λ.(T -T C ) ∀M ∈ Γ C T (M, t = 0) = T 0 (M) (2)
Where Γ P is the boundary of the polymer and Γ C the boundary of the channels. The temperature of the coolant is T C and the heat transfer coefficient, h, is related to the coolant flow rate (via Colburn correlation). Several strategies are possible to solve such problems using BEM. Mätzig [START_REF] Davis | The dual reciprocity method for heat transfer in polymer processing[END_REF] propose to use space and time Green's function. To express the domain integral in terms of equivalent boundary integrals, we introduce the DRM approximation [START_REF] Chen | Dual reciprocity method using compactly supported radial basis functions Communications in Numerical Methods in Engineering[END_REF]. The solution is defined as a series of particular solutions Tk located in each boundary nodes N n , and each internal nodes N i . We obtain Eq (3), explain in detail by Mathey [START_REF] Mathey | Automatic optimization of the cooling of injection mold based on the boundary element method Materials Processing and Design: Modeling, Simulation and Applications[END_REF]

C i T i -a Γ T.q * dΓ - Γ q.T * dΓ = Nn+N i k=1 β k C i Tik + Γ Tk .q * dΓ - Γ qk .T * dΓ (3)
Here, T and q denotes the temperature and the flux denotes, and C i is equal to 1 since the point i is inside the domain and to 0.5 on its regular boundary. The following Green's function T * and q * [1] denotes the fundamental solution of this equation. The vector β is define such as:

β = 1 a F -1 Ṫ (4)
Matrix F consists of interpolation-function values f = 1 + r at each point.

REDUCED MODELING

Usual reduced models perform the simulation of some similar problem or the desired one in a short time interval. From these solutions, the Karhunen-Loève decomposition [START_REF] Chinesta | Alleviating mesh constraints: model reduction, parallel time integration and high resolution homogenization[END_REF] can be performed, allowing to extract the most relevant functions describing the solution evolution.

The Karhunen-Loève decomposition

We assume that the evolution of a certain field T (x, t) is known. In practical applications, this field is defined at the spatial mesh nodes x i (with i ∈ {1, • • • , N}), and for some time t m = m.∆t with m ∈ {1, • • • , M}. We introduce the notation T m (x i ) for defining the vector containing the nodal degrees of freedom (temperatures) at time t m . The main idea of the Karhunen-Loève (KL) decomposition tell us how to obtain the most typical or characteristic structure φ(x) among these T m (x) ∀ M. This is equivalent to obtaining a function φ(x) maximizing α defined as:

α = P p=1 N i=1 φ(x i )T p (x i ) 2 N i=1 (φ(x i )) 2 (5)
This leads to:

P p=1 N i=1 φ(x i )T P (x i ) N i=1 φ(x i )T P (x i ) = α N i=1 φ(x i )φ P (x i ) ∀ φ (6)
where φ denotes the variation of φ(x) which can be rewritten under the form:

φ T k.Φ = α φ T .φ ∀ φ T ⇒ k.φ = α.φ (7)
We define the matrix Q containing the discrete field history, and the vector φ such that its i-component is φ(x i ). This yields to the eigenvalue problem

D = Q.Q T : Q =      T 1 1 T 1 2 • • • T 1 P T 2 1 T 2 2 • • • T P 1 . . . . . . . . . . . . T 1 N T 2 N • • • T P N      (8)
The functions defining the most characteristic structure of T P (x) are the eigenfunctions φ n (x) ≡ φ n associated with the largest eigenvalues.

A posteriori reduced model

We solve the eigenvalue problem defined by Eq (7) selecting the eigenfunctions φ n associated with the eigenvalues belonging to the interval defined by the largest eigenvalue such as Φ n 's sum is upper or equal to 99.9% of Φ N 's sum. In practice, n is much lower than N. Let us now try to use these n eigenfunctions φ n for approximating the solution. Let B be the following matrix:

B =      φ 1 (X 1 ) φ 2 (X 1 ) • • • φ n (X 1 ) φ 1 (X 2 ) φ 2 (X 2 ) • • • φ n (X 2 ) . . . . . . . . . . . . φ 1 (X N ) φ 2 (X N ) • • • φ n (X N )      (9)
We express the linear system of equations resulting from the semi-implicit thermal-model discretization as:

T m+1 = i=n i=1 ζ m+1 i φ i = B.ζ m+1 (10)

MOULD COOLING OPTIMIZATION

Each optimization iteration involves performing a BEM simulation and computing the objective and constraint functions. The optimization method allows updating the cooling channel design parameters (subject to the constraints) until a minimum of the cost function is reached [START_REF] Pirc | BEM-based cooling optimization for 3D injection molding[END_REF]. Figure 1 shows the coupling between the thermal solver and the optimization algorithm. The SQP method is designed for monoobjective optimization [START_REF] Nocedal | Numerical optimization series in operation research[END_REF]. However, pratical optimization problems almost always involves at least two objective functions. One way to proceed in such a context is to consider as cost function a weighted sum of two objectives, but this method involves choosing a weighting parameter. We rather propose here using one objective as optimization criterion, and the other as a non-linear constraint. The first criterion involves miniming the maximal temperature on the cavity surface. The second criterion aims at improving temperature uniformity. More precisely, we formulate our problem under the form:

minimize max i∈D (T i ) (11) subject to i∈D |T i -T moy | ≤ σ ( 12
)
where D is the set of discretization elements of the plastic part where the temperature T i is measured, T moy is the average of the T i 's, and σ is a user-defined temperature uniformity tolerance, fixed here equal to 4.

APPLICATION TO TWO-DIMENSIONAL CASE

It is important to note that using reduced model to optimize the cooling channel location is possible since the vector φ n (x N ) does not change when the unsteady equation source term changes, i.e. even if the position of the cooling channels changed. Optimization variables are the coordinates (X i , Y i ) of each circle center i, and the radius of them. This geometry have 7 channels, thus we have 21 optimization variables in our problem. The coolant temperature is fixed as T C = 30 • C. The reduced model permits to divide by more than 40 the CPU of each direct computation, compared to DRM. We used n = 17 in our simulations, whereas N = 288 (nodes number).

CONCLUSIONS

Our methodology uses BEM to solve the unsteady heat transfer equation during the cooling step of the injection moulding process. Simulation results are used in an optimization procedure to find the best geometry and process parameters according to a given objective function. Reduce model technique involves a Karhunen-Love decomposition leading to an optimal number of approximation functions, allowing to considerable CPU time savings (some times in the order of 40). Our preliminary test showed that our approach is viable for optimizing the design of cooling channels for injection moulding. Various objective functions can be provided by the user (either directly as a cost function or within constraints) . We presently work on more complex 3D moulds with more general parameterizations of the cooling channels.

Figure 1 :

 1 Figure 1: optimization procedure

Figure 2

 2 displays the geometry used to validate our method. Dotted circles show the initial configuration of the cooling channels, and bold circles show the optimized position.

Figure 2 :Figure 3 :

 23 Figure 2: Channels configuration before and after optimization

Table 1 :

 1 CPU time comparison

	Method	direct computation optimization
	DRM	38.2 seconds	5.6 hours
	reduction model 1 second	9.6 minutes


N Lebaal

S Puissant

F Schmidt

COMPARAISON BETWEEN TWO DESIGNS OF EXPERIMENT USING THE RESPONSE SURFACE METHOD IN ORDER TO OPTIMIZE THE VELOCITY DISTRIBUTION IN A COAT HANGER DIE

Keywords: (Optimization, DoE, Response surface method, FEM)

Balancing the distribution of flow through a die to achieve a uniform velocity distribution across the die exit is one of the most difficult tasks of extrusion die design. The objective of this paper is to obtain a homogeneous velocities distribution at the die exit. In order to keep the same geometry and avoid design and manufacture of a new die, it seems very important to control the die thermally. For this, we optimize the wall temperature of regulation of a coat hanger die in a heterogeneous way, (i.e. the wall temperature may not be constant in the entire die). The temperature of regulation of the melt enables us to locally control the viscosity, which influences the flows in the various zones.

The flow analysis results are then combined with an automatic optimisation algorithm that is based on a response surface methodology and a non linear constraint algorithm SQP with several strategies to provide a new profile of the die wall temperature distributions. Two designs of experiment are used and both results are then compared. Typically, for extrusion die design, the objective function states that the average velocity across the die exit is uniform. Constraints are used to limit and /or to control the pressure drop in the die.

INTRODUCTION

To obtain the optimal parameters during polymer extrusion processes, the numerical simulation imposed itself as an important tool by giving access to difficult orders of magnitudes of physical parameters obtain by the experiment. Moreover it improves the comprehension of the process. Far from being a substitute for the experiment, the numerical simulation allows a time saving and a substantial economy. This is obtained by reducing the number of experiments to be carried out when analyzing the process before its realization. However, obtaining the optimum parameters is often the result of a difficult work of trial and error, during which the various solutions are tested and modified.

The majority of the early works on the extrusion process use the geometrical parameters to optimise the velocity distribution at the die exit [START_REF] Nóbrega | Flow Balancing in extrusion dies for thermoplastic profiles[END_REF][START_REF] Ettinger | Parameterization and optimisation strategies for the automated design of U PVC profile extrusion dies[END_REF][START_REF] Michaeli | Computer aided optimisation of extrusion dies[END_REF][START_REF] Sun | Optimization of a flat die geometry[END_REF]. The channel geometry of a flat die should be optimized such that a uniform velocity distribution at the die exit is obtained without excessively increasing the pressure drop across the die. Within this framework Smith [START_REF] Smith | An optimisation-based approach to compute sheeting die designs for multiple operating conditions[END_REF] optimizes a flat die design to well operate at multiple temperatures. The author shows that the exit velocity distribution is influenced by the melt temperature (the material properties for the power law fluid model are varied according to the melt temperature). To simplify the optimisation approach, the lubrication approximation to model isothermal flow of power-law fluids are used. Design sensitivities needed for the gradient-based constraint optimization algorithm are evaluated with the adjoint variable method to improve the performance and accuracy of the computations.

If the non isothermal flow is not accounted while optimizing the die [START_REF] Wang | The flow distribution of molten polymers in slit dies and coat hanger die through three-dimensional finite element analysis[END_REF] the computed velocity pressure fields are expected to have large errors [START_REF] Smith | An optimisation-based approach to compute sheeting die designs for multiple operating conditions[END_REF]. The non-isothermal flow of a power-law fluid in a coat-hanger die is studied by Vergnes et al. [START_REF] Vergnes | Non-isothermal flow of a molten polymer in a coathanger die[END_REF]. The authors employed the finite difference method to solve the two dimensional equilibrium and thermal equations (flow field by considering the thermal regulation). The pressure, residence time, and temperature distribution are obtained. The authors show that the thermal regulation and temperature dependence of the viscosity have a large effect on the flow uniformity at the end of the die. Schläfli [START_REF] Schläfli | Analysis of polymer flow through coat-hanger melt distributors[END_REF] made an experimental study on the effect of the difference of the temperature between the material and the die on the velocity distribution at the exit of a Coat-hanger melt distributors. He noted that the difference in temperature has a more important effect on the final zone of the distribution channel because of the more important residence time. Within the framework, Vergnes et al. [START_REF] Vergnes | Modélisation des écoulements dans les filières d'extrusion[END_REF] were interested in the manufacture of a tube. The problem encountered during this manufacture is the presence in the section of the tube of a symmetrical extra thickness. This defect is due to the existence, in the entry of the die, of two hotter zones of melt due to the flow in the twin screw extrusion machine. These hot zones separate to flow out on each side of the punch, providing, because of their lower viscosity, more important local flows and thus, extra thicknesses at the exit. To obtain a homogeneous distribution of the flow at the die exit, instead of imposing a homogeneous temperature of regulation, they imposed a heterogeneous temperature of regulation compared to the cold and hot zones. The results show that with this heterogeneous wall temperature, they obtained a practically homogeneous distribution of thickness all around the die. This possibility of heterogeneous thermal regulation was applied and exists at the industrial level. Fradette et al. [START_REF] Fradette | Optimal-Design in Profile Extrusion Calibration[END_REF] optimize the cooling channels diameter and position for the cooling layout of vacuum calibrators used in the extrusion of PVC profiles. The required objective is to minimize the average temperature to decrease the time of cooling and to standardize the temperatures to improve quality. The methodology is based on the finite element resolution of the heat transfer conduction problem in the calibrator and the profile coupled to the Fletcher-Reeves optimization algorithm.

For complex geometries, the computational resources and time required for analysis of extrusion dies are considerable. Consequently, the use of evolution strategy algorithm to optimize the extrusion die in three-dimensional analysis is less attractive to designers because of the unacceptable turnaround time for results. The purpose of this work is to deal with this inconvenient disadvantage, by using an automatic optimisation algorithm and several strategies which permit to avoid the local minima and by this fact without increasing the computing time which becomes penalizing for three-dimensional simulations. The non isothermal flow analysis results are then combined with an automatic optimisation algorithm that is based on a response surface methodology [START_REF] Myers | Response surface methodology, Process and product optimization using designed experiments[END_REF] and a non linear constraint algorithm SQP to optimize the wall temperature profile of a coat hanger die in a heterogeneous way.

To obtain the optimal parameters during polymer extrusion processes, the numerical simulation imposed itself as an important tool by giving access to difficult orders of magnitudes of physical parameters obtain by the experiment. Moreover it improves the comprehension of the process. Far from being a substitute for the experiment, the numerical simulation allows a time saving and a substantial economy. This is obtained by reducing the number of experiments to be carried out when analyzing the process before its realization. However, obtaining the optimum parameters is often the result of a difficult work of trial and error, during which the various solutions are tested and modified.

The majority of the early works on the extrusion process use the geometrical parameters to optimise the velocity distribution at the die exit [START_REF] Nóbrega | Flow Balancing in extrusion dies for thermoplastic profiles[END_REF][START_REF] Ettinger | Parameterization and optimisation strategies for the automated design of U PVC profile extrusion dies[END_REF][START_REF] Michaeli | Computer aided optimisation of extrusion dies[END_REF][START_REF] Sun | Optimization of a flat die geometry[END_REF]. The channel geometry of a flat die should be optimized such that a uniform velocity distribution at the die exit is obtained without excessively increasing the pressure drop across the die. Within this framework Smith [START_REF] Smith | An optimisation-based approach to compute sheeting die designs for multiple operating conditions[END_REF] optimizes a flat die design to well operate at multiple temperatures. The author shows that the exit velocity distribution is influenced by the melt temperature (the material properties for the power law fluid model are varied according to the melt temperature). To simplify the optimisation approach, the lubrication approximation to model isothermal flow of power-law fluids are used. Design sensitivities needed for the gradient-based constraint optimization algorithm are evaluated with the adjoint variable method to improve the performance and accuracy of the computations.

If the non isothermal flow is not accounted while optimizing the die [START_REF] Wang | The flow distribution of molten polymers in slit dies and coat hanger die through three-dimensional finite element analysis[END_REF] the computed velocity pressure fields are expected to have large errors [START_REF] Smith | An optimisation-based approach to compute sheeting die designs for multiple operating conditions[END_REF]. The non-isothermal flow of a power-law fluid in a coat-hanger die is studied by Vergnes et al. [START_REF] Vergnes | Non-isothermal flow of a molten polymer in a coathanger die[END_REF]. The authors employed the finite difference method to solve the two dimensional equilibrium and thermal equations (flow field by considering the thermal regulation). The pressure, residence time, and temperature distribution are obtained. The authors show that the thermal regulation and temperature dependence of the viscosity have a large effect on the flow uniformity at the end of the die. Schläfli [START_REF] Schläfli | Analysis of polymer flow through coat-hanger melt distributors[END_REF] made an experimental study on the effect of the difference of the temperature between the material and the die on the velocity distribution at the exit of a Coat-hanger melt distributors. He noted that the difference in temperature has a more important effect on the final zone of the distribution channel because of the more important residence time. Within the framework, Vergnes et al. [START_REF] Vergnes | Modélisation des écoulements dans les filières d'extrusion[END_REF] were interested in the manufacture of a tube. The problem encountered during this manufacture is the presence in the section of the tube of a symmetrical extra thickness. This defect is due to the existence, in the entry of the die, of two hotter zones of melt due to the flow in the twin screw extrusion machine. These hot zones separate to flow out on each side of the punch, providing, because of their lower viscosity, more important local flows and thus, extra thicknesses at the exit. To obtain a homogeneous distribution of the flow at the die exit, instead of imposing a homogeneous temperature of regulation, they imposed a heterogeneous temperature of regulation compared to the cold and hot zones. The results show that with this heterogeneous wall temperature, they obtained a practically homogeneous distribution of thickness all around the die. This possibility of heterogeneous thermal regulation was applied and exists at the industrial level. Fradette et al. [START_REF] Fradette | Optimal-Design in Profile Extrusion Calibration[END_REF] optimize the cooling channels diameter and position for the cooling layout of vacuum calibrators used in the extrusion of PVC profiles. The required objective is to minimize the average temperature to decrease the time of cooling and to standardize the temperatures to improve quality. The methodology is based on the finite element resolution of the heat transfer conduction problem in the calibrator and the profile coupled to the Fletcher-Reeves optimization algorithm.

For complex geometries, the computational resources and time required for analysis of extrusion dies are considerable. Consequently, the use of evolution strategy algorithm to optimize the extrusion die in three-dimensional analysis is less attractive to designers because of the unacceptable turnaround time for results. The purpose of this work is to deal with this inconvenient disadvantage, by using an automatic optimisation algorithm and several strategies which permit to avoid the local minima and by this fact without increasing the computing time which becomes penalizing for three-dimensional simulations. The non isothermal flow analysis results are then combined with an automatic optimisation algorithm that is based on a response surface methodology [START_REF] Myers | Response surface methodology, Process and product optimization using designed experiments[END_REF] and a non linear constraint algorithm SQP to optimize the wall temperature profile of a coat hanger die in a heterogeneous way.

MESH SIZE SENSITIVE

In this section we examine and analyze the effect of the mesh size over CPU time and on some results: pressure in the die entrance, and the exit velocity distribution. The mesh distribution was improved by dividing the die into four zones. The first zone represents the entry of the die and the manifold. The other zones [START_REF] Ettinger | Parameterization and optimisation strategies for the automated design of U PVC profile extrusion dies[END_REF][START_REF] Michaeli | Computer aided optimisation of extrusion dies[END_REF][START_REF] Sun | Optimization of a flat die geometry[END_REF] represent the zones of low thickness of the die (land), which are respectively: the triangle, the relaxation zone, and the lips, as figure 1 indicates. On these zones we will impose various sizes of finite elements. Considering the great shape ratio (thickness over width) (fig. 2) and to respect some number of node on the thickness the anisotropic mesh is privileged (i.e. size of element which varies according to each direction). Seven mesh sizes were analyzed. In table 1, we present the characteristics of each mesh: elements size according to each direction and for the various zones, node and element numbers, CPU time and pressure drop in the die. We compare the pressure obtained by the various mesh size with that obtained by the finest mesh taken as reference. Thus, we note, according to the results obtained (table .1), that by refining the mesh and by using an adapted mesh for each zone, the difference between the pressures becomes negligible; on the other hand, the CPU time becomes more important.

In order to quantify the influence of the mesh density on the exit velocities distribution, we represent in figure .2 the results corresponding to the various mesh size. We note that the exit velocities distribution converge by refining the mesh with very close values. By analyzing more closely the curves, we remark that a stabilization of the results is established starting from mesh 4.

OPTIMISATION PROCEDURE

Optimization variables

In order to keep the same geometry and obtain a homogeneous velocity distribution at the die exit, an optimization problem is formulated, which consists in optimizing the temperature of regulation of the flat die. For this procedure, the geometry of the die as well as the mesh remains constant. On the other hand, physical parameters, which represent the wall die temperatures of regulation, will be used as optimization variables to obtain a homogeneous exit velocity distribution. 6 temperatures of regulation "T1, T2, T3, Ta, Tb and Tc" are imposed according to the width of the die (Fig. 1). Because of symmetry, we represented only half die. To decrease the number of optimization variables, we use a control curve obtained by Kriging, to allow the control of the die wall temperature profile, and to decrease the number of variables to 3.

Among the 6 optimization variables Ta, Tb and Tc are active; whereas, T1, T2, and T3 are passive. This means that T1 varies like Ta, at the distance from 0 to 200 mm of the width of the die. From the three temperatures (active variables), we can plot a curve using the Kriging interpolation to represent other die wall regulation temperatures (passive variables). The other temperatures T2 and T3 are determined according to Ta, Tb, Tc according to the curve of Kriging interpolation. The optimization problem is thus resolved with three parameters (Ta, Tb, and Tc). We can also add to it the applied constraint conditions so that the temperatures do not exceed a certain value. Initially, the temperature of regulation is constant through all the width of the die T=230°C. This constant temperature gives a non uniform velocity distribution. Velocities are more important of about 50% at the edge and weaker of about 25 % in the middle of the die compared to the average velocity. During the optimisation processes the variables permit the variation between 200°C and 280 °C. The maximum temperature of 280°C is chosen i.e. there is no polymer degradation.

Objective and constraint Function

The problem consists in determining the optimal profile of wall temperatures of this coat hanger die so that the objective function based on the exit velocity distribution is minimal; the sheet thickness will be then homogeneous. This represents a very important material profit considering the quantity of the extruded polymer sheets each year. The objective function ( ( )

x f
) represents the total relative variation on the exit velocity distribution. In order to ensure that the pressure at the entry of the die ( e P ) does not increase compared to the pressure obtained by the initial wall temperature ( e P 0 ), a constraint function ( ( )

x g ) based on the pressure at the entry of the die is used. This problem of optimization is represented by the equation ( 1).

( )

( ) ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ ≤ ≤ ≤ - = = max min 0 0 0 0 min x x x P P P x g E E x f e e e (1)
The velocity dispersion ( E ), is defined as:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = ∑ = N i i v v v N E 1 1 (2)
where 0 E and 0 P are respectively the total relative variation of the exit velocity distribution and the pressure in the initial die, N is the total number of nodes at the die exit in the middle plane; i v is the velocity at an exit node, and v is the average exit velocity which is defined as follows:

∑ = = N i i v N v 1 1
(3) min

x and max x are the lower and higher limits of optimization variables, with: The maximum temperature of 280°C is chosen so there is no polymer degradation.

Sensitive analysis

A Taguchi method [START_REF] Montgomery | Design and analysis of experiments[END_REF] is used to investigate the sensitive analysis effect of the optimization variables (regulation temperature) in the objective and constraint function. The central composite design or CCD is the most popular class of designs used for fitting a second-order model. The CCD consists of a 2k factorial with runs, 2 k axial or star runs, and central runs. The parameter, which represents the distance of the axial runs from the design centre must be specified.

The rotatability is important for the second-order-model to provide good predictions throughout the region of interest. One way to well define is to require that the model have a reasonably consistent and stable variance of the predicted response at point of interest. Because the purpose of RSM is the optimization and location of the optimum prior to running the experiment, is important to use a design that provides equal precision of estimation in all directions. For this, the CCD is made rotatable.

In this first study, in order to evaluate the effects of the optimization variables on the objective and constraint function, Taguchi design of experiments method was used to estimate the effects of the profile of the die wall temperature on the pressure drop and on the exit velocities distribution and to identify the influential parameters. This method uses a central composite design (CCD) [START_REF] Montgomery | Design and analysis of experiments[END_REF], to study the entire parameter space with only small number of experiments. In this study, three parameters were considered with five levels each. According to the Taguchi quality design concept, 15 complete simulations were chosen. The various levels are represented from 1 to 5 respectively giving the values of the minimum to the maximum of each variable. 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 "objective and constraint function " (Table 2).

Table . 2.

Response of the objective and constraint functions using CCD.

Responses

T a T b T c

E en % g Pressure [bar] 1 2 ) 1 ( α - 2 ) 1 ( α - 2 ) 1 ( α - 23.4 0.1978 5.33 2 4 ) 1 ( α + 2 ) 1 ( α - 2 ) 1 ( α - 7.04
-0.189 4.36

-0.0859 4.06 The average effects of the various variables are represented on graphs to support the discussion and to lead to the identification of those influencing to minimize the defects. Figure 4 represent the average effect of the optimization variables on the objective function (velocities distribution at the die exit). It is noted that the temperature Ta have a most important effect on the velocities distribution. This temperature is located at the medium of the die, where the velocities are lowest. For the other temperatures, Tb and Tc seem having less important effects. We notice that the combination which offers the minimum of the objective function, corresponds to the maximum of the temperature Ta, and the minimum of the temperatures Tb, and Tc. This allows to change the polymer viscosity and control the flow locally. At this place, or the air-gap is weaker, that compensates for the lack of flow on this zone. On the other hand, for the zones located on the edge of the die (more important velocities), the temperatures of regulation Tb and Tc, which are low, generate an increase in viscosity, and thus a reduction in the flow. Concerning their effects on the pressure (fig. 5) we can classify them by descending order. The greatest effect on the pressure variation is given by Ta, then Tb and finally Tc which has the weakest effect. The combination, of the optimization variables, which offers the minimum of the constraint function (pressure), corresponds to the maximum of the temperatures of regulation. That makes it possible to increase the melt temperature. Into consequent, the melt viscosity decreases, which favourite the flow with less pressure loss.

Response surface method

The objective and constraint functions are implicit, compared to the optimization parameters and their evaluation requires a large amount of computing time using three dimensional analyses. In order to decrease the evaluation number of the objective and constraint functions we use in this work the response surface method [START_REF] Myers | Response surface methodology, Process and product optimization using designed experiments[END_REF][START_REF] Lebaal | Rheological parameters identification using in-situ experimental data of a flat die extrusion[END_REF] which consists in the construction of an approximate expression of the objective and constraint functions 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 obtained using Design of Experiments technique.

Two optimization results are built (conducted) using two DoE for the construction of the R.S (figure 6) the firs one is the Enriched Composite Design (ECD) and the second is the Central Composite Design (CCD).

Kriging interpolation

Kriging is an interpolation technique [START_REF] Kaymaz | Application of Kriging method to structural reliability problems[END_REF][START_REF] Lebaal | Optimizations of coat-hanger die, using constraint optimization algorithm and Taguchi method[END_REF] which allows to model very efficiently complex functions. It is applied here to represent the response surface in an explicit form, according to the variables of optimization. The explicit relationship of the objective and constraint function can be expressed as follows:

Z(x) a (x) p (x) J T + = ~ (5) with, ( ) ( ) T m x p , , x p p(x) ] [ 1 K =
, where m denotes the number of the basis function in regression model,

T m ] a , , a [ a K 1 =
is the coefficient vector the, x is the design variables, (x) J ~ is the unknown objective or constraint interpolate function, and Z(x) is the random fluctuation. The term a (x) p T in Equation ( 5) indicates a global model of the design space, which is similar to the polynomial model in a Moving Least Squares (MLS) approximation. The second

Design of experiment

CCD ECD part in Equation ( 5) is a correction of the global model. It is used to model the deviation from a (x) p T so that the whole model interpolates response data from the function. The output responses from the function are given as:

( ) ( ) ( ) { } x f x f x f x F n K , , ) ( 2 1

=

(6) From these outputs the unknown parameters a can be estimated:

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

( ) ( ) ( ) ( ) ( ) ( )
( ) ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ - - - + ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = n 2 1 n n 1 n n 1 1 1 0 0 0 0 0 0 x , x , x , x , x x w x x w x x w x R x R x R x R R L M M M M L L L M O M L (8) ( ) 3 j i j i x x R - = (9)
A weight function of Gaussian type with a circular support is adopted for the Kriging interpolation because its derivatives with respect to the coordinates exist to any desired order. It takes the form [START_REF] Fradette | Optimal-Design in Profile Extrusion Calibration[END_REF] where

⎪ ⎩ ⎪ ⎨ ⎧ ≥ ≤ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - - - - - = w i w i w w i i r d r d c r c r c d x W χ χ * ) ) / ( exp( 1
) ) / ( exp( ) ) / ( exp( 1 ) ( 2 2 2
( ) 2
) (

1 i J J n J i x x d - = ∑ =
is the distance from a discrete node i x to a sampling point x in the domain of support with radius w r , and c is the dilation parameter.

4 w r c = is used in computation. 1 0 ≤ ≤ χ
represents the degree of the importance of the weight function. The second part in Equation ( 5) is in fact an interpolation of the residuals of the regression model a (x) p T . Thus, all response data will be exactly predicted; is given as:

( ) ( )β x r x Z T = (11)
Where ( ) ( )

( ) { } n 1 x , , , x , x R x R x r T L =
The parameters β are defined as fallow:

( )

Pa F R - = -1 β (12)

3.6.Optimization strategy

In order to reach a true global optimum (defined by the exact response) and to respect the non linear constraint function, an efficient optimization strategy using response surfaces are build:

In the first step we use the global response surface, using the DoE defined on the whole research space. The objective and constraint functions will be computed into a limited evaluation response, and we minimize it using a SQP (sequential quadratic programming method). To avoid local minima, and by this fact without increasing the computing time which becomes penalizing for three-dimensional simulations we change the initial points of the SQP method by all the points of the DoE because the evaluation of the explicit functions does not cause problems and does not take much time. When we have obtained the best minimum; successive local approximate problems based on the weight function of Gaussian type are then solved using an SQP algorithm. This weight function gives, also more importance to the responses which are closer to the minimum and less importance to the other responses, which allows makes the response surface more accurate locally. During the progression of the optimization procedure (figure 7), the region of interest moves and zooms on each optimum by reducing the search space by 1/3. In addition, We enriches then the collection { k

x } with the points by the vicinity

k x min ℜ . Giving a new collect { } { } k x k k x x min 1 ℜ ∪ = + at the stage 1 + k . At this iteration ( 1 + k ) we construct a new approximation 1 ~+ k f starting from { 1 + k x }.
The number of evaluation points considered increases because of enriching the collection { k

x } with the points by the vicinity. This carries out to an auto adaptive search space and increases the accuracy when we approach the optimum. These steps are repeated till the new profile temperature satisfies the objective and the constraints function and the iterative procedure stops when the successive points are superposed with a tolerance of ε=10 -3 .

NUMERICAL APPLICATION

The optimization strategy applied consists in using the optimization algorithm described previously, based on the response surfaces method and the Kriging interpolation. Two results are presented by using two different experimental designs:

• Case 1: Optimization is carried out by using an Enriched Composite Design of experiment (ECD).

• Case 2: We a Central Composite Design of experiment (CCD). The process of optimization is carried out and the obtimum is obtained from the 3 rd iteration; therefore a very fast convergence is obtained for the two results. These calculations take a computing CPU time of 5 days and 9 hour for the ECDs and 2 days and 9 hours for the CCD on a machine Pentium IV, 3 GHz, 1 Go RAM.

The results for the objective function (velocities distribution at the die exit) and the constraint function (pressure) are normalised (standardized) compared to the results obtained with the same geometry by imposing a constant wall temperature over all the width of the die (T w= 230°C). A first observation in table 3 shows that with the two results, we find very satisfactory solutions, and that the value of the cost function is reduced moreover than 92%. The constraints on the pressure are respected for the two solutions, as well as the velocities distribution is better. The optimal solution is obtained after three iterations for the two results. On the other hand, the computing time is more important by using the enriched composite design (ECD) (cas1). With the use of the CCD (cas2), we can obtain a reduction of CPU time of about 66%. The optimum parameters obtained are given in table 3.

Table 3. Summary of the results for the two optimization cases. The figure 8 represents the evolution of the objectives functions (E/E0) for the two optimization results (cas1 and cas2). We note that with the use of the ECD we obtain a reduction of 80% of the objective function at the first iteration, since the good interpolation obtained by Kriging. At the second and third iteration, the objective function decreases gradually until obtaining a profit of homogeneity of 93.5% compared to the exit velocities distribution obtained initially.

With the use of the CCD, the solution (total relative variation) obtained at the first iteration is higher than the initial solution obtained with a constant wall temperature. This is due to the weak evaluation points for the field of research which varies between 200 and 280°C, which makes the approximation less accurate. On the other hand, at the second iteration, the objective function decreases by 85%, then of 92.4% with the third iteration. The strategy of auto-adaptive search space makes it possible to enrich the approximation and to make it more accurate when we approach the global optimum. The figure 9 represents the relative variation of the exit velocities distribution obtained in the initial die wall temperature of regulation and after optimization for both results (using two design of experiment). Qualitatively it is clear that the velocities distribution obtained with the initial temperature (T= 230°C) gives more important relative variations, about 25% in the middle of the die and 50% at the edge. This important relative variation, result in lower velocities in the middle of the die, and more important velocities distribution on the edge. After optimization, we note that the velocities distribution is homogeneous over all the width of the die exit. The relative variations are lower than 5 % over all the width of the die. Qualitatively we find that the initial wall die temperature of regulation (Fig. 10-A) gives a not homogeneous velocities distribution. The optimal solution obtained for two optimizations using ECD and CCD (B and C), represents hotter temperatures of regulation in the middle of the die. The local flows are more important there and a compensation of the flow diminution is thus obtained compared to the initial regulation. On the other hand, the temperatures are lower on the edge to increase the polymer viscosity. This generates a reduction in the flow towards the edge. With these optimal temperatures of regulation, we could obtain a balance which makes it possible to give a homogeneous exit velocities distribution according to all the width on the die (Fig. 10 B andC). The figure 11 shows that the constraint on the pressure is respected for the two results (ECD and CCD). This constraint represents the ratio the pressure, in each iteration, and the pressure obtained with the initial temperature of regulation ( e e P P 0 / ). The optimization problem was formulated so that the ratio ( e e P P 0 / ) remains lower than 1 for each iteration, which gives a pressure decrease. We notes according to the constraints evolution according to the optimization iterations, that we obtains a profit of pressure of about 27% with the 3 rd iteration for the ECD (cas1). On the other hand, the pressure profit obtained by using a CCD (cas2) is lower, it accounts for 15% compared to the initial pressure. The pressure drop in the initial wall die temperature of regulation was specified as the limit pressure drop in the optimization algorithm. The goal was to improve the uniformity of exit velocity without increasing the pressure drop in the die. We show in the figure 12, that the pressure drop in the coat-hanger die, is lower for the two results (B, and C). The optimal temperatures of regulation obtained by using the EDC (A) are higher. This made decrease the viscosity of polymer and thus implies a weaker pressure drop (see the study of the effect of the optimization variables on the pressure figure 5). The pressure is decreased by 27%, compared to the limit pressure drop obtained by initial temperature of regulation. On the other hand, the temperatures of regulation obtained by using CCD (C) are lower compared to the result (B). This decrease the pressure drop by 16%.

Figure 13. Convergence history during an optimization process on the optimizations variables.

The graph (Fig. 13) represents the evolution of the optimization variables (temperatures of regulation) during iterations of optimization. The solutions of the two results are different; on the other hand they have the same gradient. It is noted that the wall temperatures of regulation obtained by the first optimization (cas1) are higher. The temperature in the middle of the die (T a ) is about 276°C, then it decrease according to the width of the die as illustrated on the figure 14 until the lowest temperature at the edge of the die Tc, which is equal to 223.5°C. It is noted that the difference of temperatures between «Ta, Tb» and «Tb, Tc», is respectively 19°C and 33°C. For Cas2, we notice that the temperatures are lower (fig. 15). On the other hand, the differences in temperature between the various zones are a little stronger. These differences represent 25°C and 35°C respectively between the temperatures "T a T b " and between "T b T c " (see table 3). ).This difference in temperature can generate difficulties during drawing and cooling on the outlet side of the die.

CONCLUSION

In this paper, a procedure was applied in order to optimize the distribution of temperatures of regulation of a coat hanger die. The optimisation method is based on a Surface Response Method and has a very fast convergence, which is an advantage when time-consuming flow analysis calculations are involved.

The advantage of such an approach lies in its (adjustment) facilitated of implementation to geometrical optimization. A curve of control was introduced in order to decrease the number of variables to 3 while keeping 6 temperatures of regulation. Two optimization cases were carried out according to the used experimental design. The two results show an improvement of the velocity distribution higher than 92% compared to the reference solution. The taking into account of the nonlinear constraint allowed decreasing the pressure at the end of the optimization. Nevertheless, the weakness lies in the difference of the temperatures on the outlet side of the die. This is likely to generate difficulties during the operations of drawing and cooling. In the continuation we will optimize the geometry of a cylindrical die and will carry out experimental validations.

Figure 1 .Figure 2 .

 12 Figure 1. Coat hangers die geometry.

Figure 3 .

 3 Figure 3. Wall temperatures distribution according to the width of the die (optimization variables).

Figure 4 .

 4 Figure 4. Average effect of variables levels on the objective function (exit velocity distribution).

Figure 5 .

 5 Figure 5. Average effect of variables levels on the constraint function (pressure).

Figure 6 .

 6 Figure 6. Design of experiment.

Figure 7 .

 7 Figure 7. The flowchart of the adopted optimization strategy.

Figure 8 .

 8 Figure 8. Convergence history during an optimization process on the objective functions.

Figure 9 .

 9 Figure 9. Relative variation of the exit velocities distribution with initial and optimal walldie temperature distribution using (ECD & CCD) .

Figure 10 .

 10 Figure 10. Exit velocity distributions in the initial (A), and optimized (B,C) temperature of regulation.

Figure 11 .

 11 Figure 11. Convergence history during an optimization process on the constraint functions.

Figure 12 .

 12 Figure 12. Pressure distributions in the initial (a), and optimized (b,c) temperature of regulation.

Figure 14 .

 14 Figure 14. Distribution profile of optimal regulation temperature (ECD).

Figure 15 .

 15 Figure 15. Distribution profile of optimal regulation temperature (CCD)

Figure 16 .

 16 Figure 16. Temperature distributions in the mid-plane of the initial (a), and optimized (b,c) temperature of regulation of the flat die.


R Gantois

A Cantarel

J.-N Felices

N Pirc

F Schmidt

BEM simulation of 3D updated resin front for LCM processes

The use of composite materials in large structures has increased in recent years. Aircraft industry has recently begun to investigate the field of Liquid Composite Molding (LCM) through research programs because of its ability to produce large parts at low cost. The resin motion through the fibrous perform is described using Darcy's law and mass conservation. The 3D flow is solved numerically using the Boundary Element Method (BEM). The simulations are performed using an updated lagrangian mesh in isothermal conditions. The numerical model is assessed by comparing computed results with the analytical solution in the case of radial injection from a spherical inlet.

Introduction

Composite materials have been used in the aircraft industry over the past three decades for their excellent specific mechanical properties. Large structural parts made of resin/fiber composites require specific forming processes. The traditional prepreg autoclave process produces high quality parts, but involves heavy manufacturing cost due to autoclave investment and use of high value added raw materials. As an alternative, recent technologies in Liquid Composite Molding (LCM) have become established recently under the designation of infusion processes. Liquid Resin Infusion (LRI) is one of them. It consists in forcing the liquid resin to circulate transversally into the dry preform using a distribution medium, consecutively to the vacuum application [START_REF] Correia | Analysis of the vacuum infusion moulding process: I. Analytical formulation[END_REF].

Prediction of resin flow during the infusion filling stage requires a full three-dimensional model in order to take into account the lag due to the distribution medium. Previous works revealed the adequacy of Boundary Element Method (BEM) solving the Darcy flow [START_REF] Um | A study on the Mold Filling Process in Resin Transfer Molding[END_REF][START_REF] Schmidt | Numerical Simulation of Resin Transfer Molding Using Linear Boundary Element Method[END_REF][START_REF] Soukane | Application of the level set method to the simulation of resin transfer molding[END_REF]. It was performed assuming a twodimensional flow, neglecting the impregnation through the thickness.

The present paper is focused on a three-dimensional numerical modeling of the Newtonian isothermal flow using constant Boundary Elements. Combining Darcy's law with incompressibility equation leads to Laplace equation which is solved using BEM at each step time on an updated lagrangian mesh. Preliminary results are presented, i.e. BEM simulations performed in the case of central injection. Computed results are compared with an analytical solution.

Governing equations

Modeling assumptions

The liquid resin flow through homogenous and not deformable porous medium (i.e. fibrous reinforcement) is modeled by the phenomenological Darcy's law [START_REF] Darcy | Les Fontaines Publiques de la Ville de Dijon. V. Dalmont[END_REF] together with incompressibility equation. The filling stage is assumed to be isothermal.

Darcy's law

The average velocity v r is proportional to the pressure gradient according to the following equation:

[ ] p K v ∇ - = r r µ (1)
where [ ] K is the permeability tensor, p the resin pressure and µ the liquid resin viscosity. Combining Equation 1 and incompressibility equation gives:

[ ] 0 . =         ∇ - ∇ p K r r µ (2)
For sake of simplicity, we assume in this paper that the principal orthogonal axis of the reinforcement coincides with the reference coordinates system, such as [ ] K reduces to diagonal matrix. Equation 2 can then be reduced to Laplace equation in the isotropic equivalent domain by stretching the coordinates. In the Cartesian coordinate system, the equivalent coordinates e x r may be deduced from the physical coordinates x using the following transformation:

i i e e x K K x i = (3)
where i K is a principal value of [ ] K and e K the equivalent permeability given as

3 3 2 1 . . K K K K e =
. In the isotropic equivalent domain, Equation 1 and Equation 2 reduce simply to:

p K v e e ∇ - = r r µ (4) 0 = ∆p (5)

Analytical solution

Following the approach introduced by Adams et al. [START_REF] Adams | Radial penetration of viscous liquid into a planar anisotropic porous medium[END_REF], an analytical solution is used to calculate the evolution of the flow front versus time. Computations are performed in the isotropic equivalent system, using the equivalent permeability e K and the spherical injection gate of radius e r 0 . A constant pressure injection 0 p is applied, so that the macroscopically resin front shape is spherical of radius ( )

        - - = e e e f e f f e f r r r p p K v 1 1 0 2 0 µε (7)
The evolution of the flow front is governed by the following differential equation: ε µ

( )         - - = =
- = +         -         (9)
Equation 9 gives immediately the flow front which is spherical in shape for an isotropic reinforcement. For the orthotropic case, the reciprocal transformation is applied to give the three semi-axis of the flow front which is ellipsoidal in shape.

Numerical simulation

Boundary Element Method

The calculation domain is delimited by its boundary

U q p Γ Γ = Γ
. A known value of pressure p (resp. normal pressure gradient q ) is prescribed on the boundary surfaces p Γ (resp. q Γ ). Equation 4is multiplied by the Green function * p . Integrated by parts twice over the calculation domain and using

Green's theorem leads to the well-known Somigliana equation [START_REF] Brebbia | Boundary elements: an introductory course[END_REF][START_REF] Pirc | BEM-based cooling optimization for 3D injection molding[END_REF]:

π θ 2 with * * = Γ = Γ + ∫ ∫ Γ Γ i i i c d q p d pq p c (10)
where i p is the value of the pressure at a point i M located on the boundary (which is supposed to be smooth), * q the pressure gradient associated with * p and θ is the internal angle of the corner in radians.

For a three-dimensional domain, * p and * q are given as: Introducing the N*N matrix H and G , Equation 12 is then transformed into:

ij 1 1 G q H p N j j ij N j j ∑ ∑ = = = (13) ∫ ∫ Γ Γ Γ = Γ + = j j d p G d q H ij ij
and 2 1 where * * ij δ Finally, Equation 13 is reordered to take the form of a linear system as:

F AX = (14)
where X is the vector of N unknowns, A is a N*N matrix obtained by reordering H and G . F is the know vector computed from the boundary conditions, H and G matrix.

Fluid domain updating

The fluid domain is continuously expanding as it flows, constraining to update the resin flow front. The filling stage is regarded as a succession of quasi-steady states, so that the geometry is redefined at each time step. Unknown pressures p and normal pressure gradients q are calculated using Equation 14. Nodes on the moving boundary are relocated using a backward Euler integration scheme:

n n v t t x t t x ) . ( ) ( ) ( ∆ + = ∆ + (15)
Combining Equation 15and Equation 4 yields to:

) ( ) ( ) ( n q K t t x t t x e µ ∆ - = ∆ + (16)
The procedure updates the resin flow front, but induces confined mesh degeneration. To overcome this difficulty, an additional adjustment is performed relocating the nodes in contact with the mold.

Algorithm

The implemented algorithm is summarized below:

• Initial computation domain is meshed.

• Stationary boundary conditions p or q are applied on each BEM element.

• Nodes coordinates are transformed according to Equation 3.

• Normal pressure gradients q are calculated according to Equation 14.

• Normal pressure gradients are prescribed on each node averaging values from surrounding elements.

• Nodes on the moving boundary are relocated according to Equation 16. Nodes in contact with the upper mold wall are relocated according to a homogenization procedure. • New normal pressures p and normal pressure gradients q are calculated next steps. The three last steps are repeated until the allotted mold-filling time is reached. • Nodes coordinates are transformed according to the reciprocal transformation.

Validation of numerical simulation 4.1 Boundary conditions

The simulation consists in propagating the resin front into an anisotropic reinforcement characterized by its permeability tensor [ ] K and its porosity ε . A constant pressure 0 p is prescribed at the spherical inlet of radius 0 r . On the flow front, the pressure f p is constant and set by the vacuum pump. Processing and rheological parameters assigned are referenced in Table 1. They keep close to the LRI process conditions but the distribution medium is not taken into account as a first approach.

Figure 1 shows the boundary conditions prescribed on the boundary of the initial fluid domain. The non penetration condition is assumed assigning 0 = q at the upper mold wall. Figure 1: Boundary conditions Table 1: Rheological and processing parameters

Comparison with analytical solution

For symmetry reasons, the simulation is performed on a quarter domain meshed using 540 elements as shown in Figure 2. The step time is set at 10 -2 s. CPU time for the complete simulation is around one hour on Intel® Core™2 Duo CPU using Matlab® environment. The computed flow fronts are in fair agreement with the analytical one as shown in Figure 3. At 10 s elapsed time, the deviations between the numerical and analytical flow fronts positions do not exceed 1%. Prolongating the simulation until 60 s gives deviations around 3%, due to the gradual mesh degeneration as the iterations go along. Figure 4 illustrates the flow front evolution during the filling stage, comparing the front shapes at t=0.5 s and t=10s.

The initial mesh generation is performed on the equivalent isotropic domain. In order to assess the mesh effect on the numerical method, some simulations have been carried out meshing the physical domain. Deviations are all the more important that the permeability values are different. This is due to the transformation from Equation 3 which stretches the mesh, degenerating elements before any calculation.

Conclusions

The resin flow front propagation is accurately simulated using a backward time-stepping algorithm. This technique is efficient for short distances of impregnation, giving a high accuracy of the front position. However, it is unsuitable for larger distances, since the calculated position becomes inaccurate as soon as elements are degenerated. Including remeshing procedures is possible but could be time consuming.

Finally, the implementation of the present numerical method becomes tedious for dividing and merging resin fronts resulting from interior obstacles or multi-injection ports. Future works will investigate the combination of Boundary Element Method with fixed grid techniques such as the Level Set method.

  such configuration the resin pressure varies only in the radial direction according to: the expression of the gradient pressure into Darcy's law and introducing the system porosity ε gives the radial component of the front velocity e f v :

  time dependant flow front position can be deduced integrating the previous equation during the filling stage. This leads to the following solution:

  r is the distance from the point i M of application of the Dirac delta function to any point under consideration. Meshing the boundary into N constant boundary elements and applying Equation 10 gives:

K

  (10 -11 m²) µ (Pa.s)

Figure 2 :Figure 4 :

 24 Figure 2: Initial meshed domain Figure 3: Comparison with analytical solution

7-10 September 2009 -Tarbes -France ICSAAM 2009

Acknowledgments

This work was supported by DAHER-Socata and CRCC-Resource and Competencies Center for Composite Materials. The authors are especially grateful to M. Niquet for its assistance during this works.


F M Schmidt
email: schmidt@mines-albi.fr

S Nakouzi

J Pancrace

Y Le Maoult

F Berthet

CURING SIMULATION OF COMPOSITES COUPLED WITH INFRARED HEATING

Keywords: polymer composite, infrared heating, degree of cure, epoxy resin, graphite/epoxy, modeling

Because of higher specific strength and stiffness, low weight, and good resistance to corrosion, the use of composite materials in aerospace structures has increased. Aircraft industry has recently begun to investigate Liquid Composites Molding techniques (LCM) through research programs because of its ability to produce large parts at a low cost. In this paper, we have not addressed the filling step during which the resin flows through fibrous media, but we investigate the numerical simulation of curing reinforced RTM-6 by infrared heating. A ray tracing in-lab software called RAY-HEAT is used to compute the radiation energy due to the interaction between reinforced RTM-6 and lamps. Finite element based program COMSOL Multiphysics™ has been used to simulate the curing process. Thermochemical model has been implemented in order to compute reaction rate as a function of reaction temperature and degree of conversion using a cure kinetic model .

INTRODUCTION

Using thermal radiation to cure the composite is more efficient, effective, and clean than standard thermal energy sources. Thermal radiation was successfully used to produce a high quality polymer matrix in less time [START_REF] Chern | On-Line Processing of Unidirectional Fiber Composites Using Radiative Heating: I. Model[END_REF]. The infrared radiation was studied and coupled with the energy equation and the exothermic heat released during the curing process in order to predict the composite temperature as a function of time and degree of cure.

INFRARED INTERACTION

Because the graphite fibers embedded in the epoxy resin are very strong infrared absorbers over broad spectral ranges at all wavelengths, the composite has strong absorption bands in some portions of the spectrum; we have neglected epoxy absorption in this paper [START_REF] Chern | On-Line Processing of Unidirectional Fiber Composites Using Radiative Heating: I. Model[END_REF][START_REF] Chern | On-Line Processing of Unidirectional Fiber Composites Using Radiative Heating: II. Radiative Properties, Experimental Validation and Process Parameter Selection[END_REF][START_REF] Kim | Transient Thermal Modeling of In-Situ Curing During Tape Winding of Composite Cylinders[END_REF][START_REF] Chern | Thermal Analysis of In-Situ Curing for Thermoset, Hoop-Wound Structures Using Infrared Heating: Part II -Dependent Scattering Effect[END_REF].

BOUNDARY CONDITIONS

The composite is assumed to be a block with a thickness range of 1-25mm. It is composed of the superposition of a sequence of carbon plies embedded in an epoxy matrix [START_REF] Cunningham | Prediction of the temperature profile with a composit sheets during pre-heating[END_REF].

Figure 1:

.Schematic presentation of the composite [START_REF] Zebdi | An Inverse Approach Based on Laminate Theory to Calculate the Mechanical Properties of Braided Composites[END_REF].

The infrared radiation is absorbed by the carbon ply present in the top surface of the composite. The boundary conditions taken in COMSOL Multiphysics TM are a surface incident heat flux on the top surface. There is also natural heat convection between this surface and the air present in the oven. The other boundaries are assumed to be adiabatic: the bottom surface is in contact with a mould and there is no heat exchange between the mould and the composite. The other sides of the composite are so thin in comparison with the top that they can be considered adiabatic.

ENERGY EQUATION

The composite transient temperature can be evaluated by the Fourier's modal equation assuming an instantaneous equilibrium temperature between the resin and the fibers at each time: [START_REF] Bailleul | Identification des propriétés thermiques de composites fibres de verre / résines thermodurcissables Application a l'optimisation des procédés de moulage[END_REF][START_REF] Ruiz | Thermomechanical Properties during Cure of Glass-Polyester RTM Composites: Elastic and Viscoelastic Modeling[END_REF] (1) Where = ultimate heat of reaction, = degree of cure, and the , , thermal conductivity of the composite are expressed using the rule of mixture: [START_REF] Ruiz | Thermomechanical Properties during Cure of Glass-Polyester RTM Composites: Elastic and Viscoelastic Modeling[END_REF][START_REF] Ruiz | Multi-criteria thermal optimization in liquid composite molding to reduce processing stresses and cycle time[END_REF][START_REF] Balvers | Determination of Cure Dependent Properties for Curing Simulation of Thick-Walled Composites[END_REF] (2)

(3) (4) (5) (6)

Where

= transverse thermal conductivity of the fibers, weight fraction, = volume fraction, the subscript stands for the fibers, and for the matrix. The system of equations described above can be solved at each time step in order to calculate the distributions of temperature and degree of cure and their evolution with time.

KINETIC MODELING

GLASS TRANSITION TEMPERATURE

The glass transition temperature is modelled using the DiBenedetto equation: [START_REF] Ruiz | Multi-criteria thermal optimization in liquid composite molding to reduce processing stresses and cycle time[END_REF] Where and = glass transition temperature of unreacted and full-reacted resin, = material constant. These values were deduced from the literature [START_REF] Ruiz | Thermomechanical Properties during Cure of Glass-Polyester RTM Composites: Elastic and Viscoelastic Modeling[END_REF].

KINETIC MODEL

Kamal and Sourour [START_REF] Bailleul | Identification des propriétés thermiques de composites fibres de verre / résines thermodurcissables Application a l'optimisation des procédés de moulage[END_REF] have shown that the following model adequately describes the cure kinetics of an epoxy resin:

Where and = rate constants with an Arrhenius type of dependence with temperature and and = catalytic constants. Although the model of Kamal and Sourour [START_REF] Zebdi | An Inverse Approach Based on Laminate Theory to Calculate the Mechanical Properties of Braided Composites[END_REF] contains several essential features, it ignores the influence of the glass transition temperature on the diffusion of the macromolecules.

This model [START_REF] Bailleul | Identification des propriétés thermiques de composites fibres de verre / résines thermodurcissables Application a l'optimisation des procédés de moulage[END_REF] has been extended in order to reflect the effect of vitrification in the kinetic model of resin transfer molding. This model was also used by [START_REF] Ruiz | Thermomechanical Properties during Cure of Glass-Polyester RTM Composites: Elastic and Viscoelastic Modeling[END_REF] [11] Where = diffusion factor, = diffusion rate constant, = constant having the order of one, = glass transition temperature, the free volume fraction of the polymer [START_REF] Nguyen | Systèmes époxy-amines incluant un catalyseur externe phénolique : Cinétique de réticulation-vieillissement hydrolytique[END_REF],

= the free volume fraction for a temperature equal to , = the thermal expansion coefficient for a temperature above . For an amorphous polymer [START_REF] Balvers | Determination of Cure Dependent Properties for Curing Simulation of Thick-Walled Composites[END_REF] which is the case of our composite, and .

KINETIC PARAMETERS

Kinetic parameters of RTM6 resin are calculated using linear regression algorithms by fitting model equations with the Differential Scanning Calorimeter experimental data. These experimental DSC data were scanned isothermally for different temperatures: 140, 160, 170 and 180°C. In the DSC, the rate of heat generation is recorded with time. As is constant for a given curing process, [START_REF] Bailleul | Identification des propriétés thermiques de composites fibres de verre / résines thermodurcissables Application a l'optimisation des procédés de moulage[END_REF][START_REF] Ruiz | Thermomechanical Properties during Cure of Glass-Polyester RTM Composites: Elastic and Viscoelastic Modeling[END_REF][START_REF] Elsawi | Rheological and kinetic behavior of double-walled carbon nanotubes filled epoxy resin[END_REF] the rate of the degree of cure is calculated by:

The degree of cure is obtained by integrating the DSC data and normalizing with .

FITTED KINETIC PARAMETERS

The figure 2 The estimated kinetics parameters from the fitting curves used in our simulations are listed in table 1 below. The carbon fibers embedded in the epoxy are very strong absorbers in that they absorb infrared energy and convert this energy to heat the epoxy matrix surrounding the carbon fibers. Thermal radiative boundary conditions are imposed on the top surface of the composite and the curing temperature is controlled through the infrared emitters in order to optimize a uniform surface heat flux. The emitters used in our project are halogen lamps given from TOSCHIBA LIGNTNING. Until now, the composite properties were taken from the literature [START_REF] Chern | On-Line Processing of Unidirectional Fiber Composites Using Radiative Heating: I. Model[END_REF][START_REF] Ruiz | Multi-criteria thermal optimization in liquid composite molding to reduce processing stresses and cycle time[END_REF]; others are assumed to be constants in order to run some initial simulations of a curing process of RTM6 as a matrix combined with a fibrous reinforcement in the form of carbon plies.

In this numerical example, a high heat flux of was provided for the first 6 minutes in order to initiate and accelerate resin curing, and for the rest of time, the surface heat flux taken was These heat flux values was imposed in a way to cure the composite in a less time and do not exceed 220°C in temperature (the degradation temperature of the matrix).

NUMERICAL RESULTS

Using COMSOL Multiphysics, the energy equation was coupled with the exothermic reaction in order to predict the temperature distribution in the composite and the degree of cure. The figures 4 and 5 below show the temperature and the degree of curing respectively evolution with time at the center of the composite. These numerical simulations correspond for a 1mm thickness composite; there is no important variation of the temperature or degree of cure with the thickness. From these simulations, we can deduce that the cure of the composite is accomplished in a less time than the others modes of heating. In the future work, we need to validate our numerical simulations with experimental infrared curing.

CONCLUSIONS

The cure model of resin is investigated in this study. A modification of the Kamal-Sourour cure model is given to account for the diffusion effect. Kinetic parameters of the cure model based on DSC Experiments were estimated using optimization algorithms. A thermalkinetic model has been developed in COMSOL Multiphysics to simulate the cure process of pure resin with assumed material properties. In-lab software developed to predict the infrared heat flux on the top surface of the composite will be coupled with COMSOL in order to simulate the whole infrared curing process of the composite.

  below shows the temperature dependence of factors and which exhibits an Arrhenius type of temperature dependence.

Figure 2 :

 2 Figure 2: Reaction rate constants k1 and k2.

Figure 3 :

 3 Figure 3: Schematic presentation of the composite face to the infrared emitters.

Figure 4 :

 4 Figure 4: Temperature evolution with time at the center of the composite

Figure 5 :

 5 Figure 5: Degree of cure evolution with time at the center of the composite

Figure 6 :

 6 Figure 6: Temperature evolution with thickness and with time

Figure 7 :

 7 Figure 7: Degree of cure evolution with thickness and with time

Table 1 :

 1 Estimated parameters

	parameter
	[	]	4.5E6
	[	]	74.69E4
	[	]	1.3E6
	[	]	58.37E4
	[	]	85.2E28
	[	]	2.16E5
	[-]		1.23
	[-]		1.36
	[-]		0.227
	[J/kg]		4.30E5
	R [J/(mol-K)]	8.314

ACKNOWLEDGEMENT

The authors would like to thank Mr. Goyot D. for his support. Grateful acknowledgement to TOSCHIBA LIGHTNING for their financial support and their attention to this work.


Sawsane Nakouzi
email: sawsane.nakouzi@mines-albi.fr

Johann Pancrace

Fabrice Schmidt

Yannick Le Maoult

Florentin Berthet

Simulation de la polymérisation par infrarouge du composite

Résumé -Une méthode récente d'élaboration des composites est la polymérisation par infrarouge. Dans cette technique, le renfort fibreux, est imprégné de résine à l'aide de la technique d'infusion de résine sous vide (en anglais Liquid resin infusion). Ce composite est placé dans un four IR, face aux lampes halogènes. Dans cette étude, les propriétés radiatives des systèmes carbone/époxy ont été déterminées expérimentalement à l'aide d'un spectromètre IR. Une première approche a été effectuée également sur le système verre/époxy. Les mesures expérimentales ont montré que le transport radiatif dans un renfort de carbone est surfacique: la pénétration des rayons infrarouge dans l'épaisseur est négligeable, d'où chaque rayon IR est absorbé par le premier pli de carbone présent à la surface du composite. Au contrario, le transfert radiatif est volumique pour le système verre/époxy. Dans cet article, on présente les simulations du processus de cuisson par infrarouge du système verre/époxy : Un logiciel développé au laboratoire nommé RAYHEAT basé sur la méthode de lancer de rayons nous permet de prédire la distribution du terme source radiatif lors de l'interaction IR avec un matériau semi transparent. Ce logiciel a été modifié afin de prendre en compte la divergence du terme source radiatif due à l'anisotropie des fibres de verre. Ce terme source est exporté vers le logiciel commercial COMSOLMULTIPHYSICSTM où l'équation de chaleur est couplée avec la cinétique de réaction de la résine. Dans des résultats présentés on a simulé la cuisson du système carbone/époxy où le transport radiatif est surfacique et nous proposerons des résultats préliminaires sur le chauffage du système verre/époxy.

Nomenclature

Introduction

L'amélioration des performances des composites à hautes performances passe par de nouvelles méthodes de mise en oeuvre de ces matériaux. L'utilisation de l'infrarouge [1] pour la polymérisation des résines époxy couramment utilisées dans les applications aéronautiques, est une méthode de mise en oeuvre qui nous permet d'assurer les performances mécaniques finales élevées du composite, dû à la flexibilité de ce mode de chauffage. Ce mode s'applique aussi au préchauffage du renfort [5]. Par contre, la prédiction de la distribution de température au cours de cuisson par infrarouge, est très importante, la résine ayant une température de dégradation thermique. D'où la nécessité d'un logiciel de calcul nous permettant de prédire cette température en fonction du temps. Le laboratoire à l'ICA Albi a fait plusieurs recherches [7] destinées a une meilleurs compréhension de l'interaction des rayonnements IR avec les matériaux. Un logiciel nommé RAYHEAT basé sur la méthode de lancer de rayons a été développé. Ce logiciel permet de prédire la distribution du terme source radiatif dans un matériau semi-transparent et du flux radiatif surfacique dans le cas d'un matériau opaque. La validation de ce logiciel par une solution analytique de facteur de forme [START_REF] Leuenberger | ASME Annual Meeting[END_REF] Figure 3 Présentation schématique du four infrarouge

Lancer de rayons

Le principe général de la méthode de lancer de rayons, consiste a simuler l'interaction des rayons émis par les lampes infrarouges et chaque élément présent dans le four (réflecteur, composite). Cette méthode est basée sur un calcul purement géométrique [7]. Pour chaque rayon, on choisi un point P sur la surface extérieure du filament tungstène (ayant une forme cylindrique) de la lampe IR. Ensuite, la direction de ce rayon est aussi choisie aléatoirement, dans le domaine de l'hémisphère formé par la normale sortante du point P. La figure 4 présente schématiquement cette méthode.

Interaction rayon du rayon IR avec le composite

Des mesures des propriétés radiatives ont été réalisées afin de bien comprendre l'interaction du rayon avec le composite dans le four. Des mesures en transmission pour les couches recouvrant le composite (figure 5,6) ont été faites pour déterminer la fraction de rayons IR absorbée par le composite. D 'autre part, le renfort en carbone est parfaitement opaque au rayonnement infrarouge, sur la figure 7 on observe les mesures en transmission pour un tissu de renfort en carbone d'épaisseur

. Ces mesures montrent que le transport radiatif à travers le renfort est surfacique, d'où les rayons IR sont absorbés par la surface du composite. D'où l'intérêt de mesurer la réflexion de la séquence des couche de film de recouvrement, tissu d'arrachage et renfort fibreux (figure 8).

Conditions aux limites

Le composite étant considéré comme opaque dans l'infrarouge, le rayonnement émis par les lampes est absorbé par le premier pli de carbone présent sur la surface supérieure du composite.

Par conséquent, la distribution surfacique du flux radiatif générée par RAY-HEAT est prise comme condition thermique aux limites sur la face supérieure du composite dans COMSOL MULTIPHYSICS. En même temps cette face est soumise à un transfert par convection naturelle dû au contact avec l'air présent dans le four. La face inférieure avec une plaque en téflon, celle ci est à son tour en contacte avec une plaque aluminium constituant le moule. La température du moule est une condition thermique aux limites. Les autres faces sont considérées adiabatiques ayant isolé thermiquement les bords du composite.

Equation de chaleur

La prédiction spatiale et temporelle de la température et du degré de polymérisation nécessite la résolution de l'équation de la chaleur, en supposant un équilibre thermique entre le renfort fibreux et la matrice à chaque instant [START_REF] Nguyen | Systèmes époxy-amine incluant un catalyseur externe phénolyque: cinétique de réticulation-vieillissement hydrolytique[END_REF]. La répartition de notre flux radiatif n'étant pas toujours homogène, elle est liée directement à la configuration géométrique du four, d'où la résolution tridimensionnelle de l'équation de la chaleur :

(1) Cette équation est fortement couplée avec l'équation de la cinétique. Ce couplage dans COMSOLMULTIPHYSICS nous permet de prédire la distribution de température ainsi que le degré de réticulation dans le composite au cours de la polymérisation [START_REF] Nguyen | Systèmes époxy-amine incluant un catalyseur externe phénolyque: cinétique de réticulation-vieillissement hydrolytique[END_REF].

(2)

(3) Les propriétés cinétiques et thermo physiques du composite ont été prises dans la littérature [START_REF] Ruiz | [END_REF]et déduites de la loi des mélanges.

Mesures expérimentales

Le montage expérimental consiste a mesurer la température du composite au cours de la polymérisation. Des thermocouples type k, de diamètre 80µm, sont placés dans l'épaisseur du composite à travers les couches de renfort.

Conclusion

Une nouvelle méthode de mise en oeuvre des composite c'est la polymérisation dans un four infrarouge a été présentée dans ce papier. Le logiciel développé au laboratoire sous MATLAB, basé sur la méthode de lancer de rayons a été utilisé afin de prédire la distribution du flux radiatif à la surface du composite. Ce flux, exporté vers COMSOL a été couplé avec l'équation de chaleur et l'équation de cinétique du composite afin de déterminer la distribution de la température dans le composite. Des mesures des propriétés radiatives du composite ont été présentées. Des mesures expérimentales nous nt permis de valider le processus de cuisson du composite. Dans le futur, nous allons simuler le chauffage par infrarouge sur une géométrie industrielle.

Figure 1 Figure 2

 12 Figure 1 Four infrarouge

Figure 4 Figure 5 Figure 7 Figure 8

 4578 Figure 4 Présentation schématique de la méthode de lancer de rayons

Figure 9

 9 Figure 9 Conditions thermiques aux limites

Figure 70

 70 Figure 70 Positions des thermocouples selon l'épaisseur du composite La figure 11 présente le flux radiatif reçu par la face supérieure du composite en son centre.La figure12présente la température de l'air en fonction du temps. Le coefficient de convection naturelle de la figure13a été calculé en fonction de la température de la face supérieure du composite, de la température du four et de la longueur du composite, en utilisant une équation semi empirique dans la librairie COMSOL.

Figure 8 Figure 10 6 )Figure 12 Figure 14

 81061214 Figure 8 Flux radiatif au centre du composite Figure 9 Température de l'air dans le four

Remerciements

Nous tenons à remercier la société TOSHIBA LIGHTING pour leur support et l'intérêt qu'ils apportent à ce travail. Nous remercions aussi B. Cosson et M. Bordival pour leur contribution au développement du logiciel RAY-HEAT.