Michel Bellet

Marie Hélène Vantal

Dr Bernard Monasse

Thermal effects in the numerical simulation of the thermoforming of multilayered polymer sheets

This paper mainly treats the thermal effects during the thermoforming process while most of the previous analyses consider an isothermal deformation. A non isothermal three dimensional finite element model of the thermoforming process is proposed. It couples the thermal equations in the thickness and mechanical equations on the mean surface of the sheet.

The mechanical resolution is done by a finite element method using a membrane approximation. The deformation is driven by a pressure difference through the sheet. The thermal resolution uses a one dimension finite element method in the thickness with convection or conduction at the surface and dissipation of mechanical energy. The polymer cooling is very efficient during the contact with tools. The coupling is done by the thermal dependent rheology. The respective contribution of friction and thermal effects in the thickness of the part during the process are discussed. The model also considers a possible multilayered material, with specific rheological parameters inside each layer. The rheology of a polystyrene was measured under elongation as a function of temperature, strain and strain rate and described by a visco plastic law. The predictions of the model were compared with measurements on an instrumented thermoforming machine and with the local thickness of axisymmetrical parts and with 3 D parts thermoformed with the same polystyrene.

INTRODUCTION

The thermoforming process consists of heating a polymer sheet and shaping it inside a mold ([1], see fig. 1). The deformation results from an applied differential pressure between the upper and lower side of the sheet, varying with time, and eventually coupled with the use of a moving plug. The forming of thin products made of various polymer materials is carried out by this method : rubbery amorphous, solid semi crystalline or multi layer composites. The heating temperature depends upon the selected polymer. At low temperature, the forming is limited by the high rigidity of the sheet. On the contrary, at high temperature, the sheet deforms by gravity and the forming is very difficult to control. Hence, thermoforming is performed above the glass transition temperature (T g ), in the rubber like behavior domain, for amorphous polymers, and close to the melting temperature for semi crystalline polymers.

One of the main problems of the process is the thinning in the corners of the parts, which leads to a decrease of the mechanical properties of the shaped components. The optimization of the final thickness profile is generally achieved by trial and error, changing the design of the component, the polymer material, and the process parameters such as the heating temperature distribution, the mold temperature, the pressure curve, or using a plug for deep parts. The numerical modeling of the shaping stage should then result in a more efficient optimization of the process : this is the point on which the present paper is focused.

Various numerical simulation models have been proposed [2 7]. Concerning the application to three dimensional forming, they are generally based on the membrane mechanical approximation, associated with the finite element method. Those computations use either hyperelastic [3, 5 7] or viscoelastic [2,4] constitutive equations, without any heat transfer coupling, and a sticking contact hypothesis at the polymer tool interfaces. It should be pointed out that this latter assumption has no experimental support. Most of the time, it is used just because the computation is isothermal and cannot account for the decrease of the polymer temperature and its "freezing" after tool contact. These models are unable to predict the effect of the mold temperature. Another thermo mechanical coupling is the self heating source term due to the high strain rates As the thermal diffusivity is low and the behavior highly temperature dependent, it seems essential to couple the heat transfer and mechanical models.

In the present paper, on line measurements on a thermoforming machine are done and compared with the kinematics predicted by the three dimensional finite element code. Local thickness measurements on the final product check another prediction of the thermo mechanical model. Concerning the contact conditions with plug and mold, either sliding contact with Coulomb's friction law or sticking contact can be considered. The paper is focused on the thermal evolution effects during the forming, and on its role in the process compared to mechanical effects such as contact conditions. Then a first approach to multilayered sheet forming is proposed.

MATERIAL AND EXPERIMENTAL METHODS

Thermoforming device

The experiments were done on two different Kiefel thermoforming machines : an instrumented machine with an axi symmetrical mold, and an industrial in line machine including multi imprints with non symmetrical shapes.

On the industrial Kiefel machine the infra red heating and thermoforming sequence is fixed. The forming is driven by pressure without any operating plug. This industrial machine is not instrumented, the maximum pressure applied is known (P = 0.6 MPa) but not the pressure variation during the forming process. Thus it is not possible to check the kinematics of the film during the process. It is only possible to measure the final profile along various directions of a complex geometry. With this device the robustness of the three dimensional finite element code and its ability to predict a complex shape can be tested. A polypropylene is thermoformed in this machine and the rheology under uniaxial tension was measured previously [8].

The pilot machine Kiefel KD20/25 includes two positions : a heating zone and a forming position. The different steps of the sequence heating transfer forming are controlled separately. The heating of the film is done by two plates including twelve quartz infra red heaters each. The thermal homogeneity of the film is controlled by the individual adjustment of the power of the heaters. The temperature is measured in the center of the film by an infra red sensor up to the final temperature, then the film is automatically transferred to the forming zone. The forming can be done by a differential pressure with pressure on the top and air pumping in the mold (Fig. 2). The aluminum mold can be heated. It is composed of two parts : a basically truncated conical cavity and an axial removable insert. The mold geometry is mainly characterized by the depth/radius ratio d/r = 60 mm / 70 mm to produce a significant deformation of the film. The insert must anticipate the film contact and then the contact contribution on further deformation. The mold is instrumented with pressure, temperature and contact gauges in various locations of the mold. The temperature is measured on the film with an infra red sensor fixed on the mold just after the transfer from heating zone to take into account possible cooling. This temperature is the initial thermal condition for the forming process. The pressure is measured on the top of the cover and in the outer corner of the mold (Fig. 2). The contact of the polymer with the mold is detected by four sensors located in critical regions including the first contact, on the insert, and the last ones, in the corners of the mold. The pressure sensors were calibrated in static condition at P = 0.1 MPa and a maximum error of δP = 5 10 4 MPa was measured on each sensor. Then it is possible to measure precisely (δP = 10 3 MPa) the pressure difference through the sheet during the thermoforming process.

Figure 3 presents a typical measurement of pressures in the cover P o and in the mold P 1 , respectively, and the contact time of the various contact gauges (C1, C2, C3). The differential pressure ∆P which acts on the deformation of the film is deduced from pressure measurements, ∆P = P o P 1 (Fig. 3). The measured differential pressure ∆P or cover pressure P o will be applied in the thermo mechanical model. The time contact measured and predicted by the model will be compared.

The forming parameters are as follows:

T m = 20°C or 60°C for the mold, T air = 20°C, d(P o )/dt = 0.125 ; 0.36 or 0.6 MPa.s 1 up to 0.1, 0.2 or 0.4 MPa, respectively, ∆P = f(t) are deduced from on line measurements, T s = 110, 130 or 140°C for the sheet.

Polymer and elongational rheology

An atactic polystyrene provided by the ELF ATOCHEM Company was mainly used during these experiments. It is a blend of 50% impact polystyrene with 50% crystal polystyrene, M w = 1.9 10 6 g/mol ( M w / M n = 3) and M w = 2.2 10 6 g/mol ( M w / M n = 2.6) respectively. The polymer was provided as extruded polymer sheets 0.3 mm thick and 20 cm large.

The rheology under uniaxial elongation of polystyrene is measured using a sheet in an isothermal condition under a constant strain rate ε . The experiments were done with a hydraulic tensile machine, Instron 1341, in an oven with a transparent door to observe the sample as it is heated by convected air. A sample 20 cm long and 15 cm wide is cut in an extruded cast film with respect to the processing conditions. It is marked with a regular grid ∆L o = 1 mm and ∆l o = 5 cm, respectively along the length and the width directions. Its distortion during the experiment gives the local Hencky deformation ε = ln(∆L / ∆L o ). The sample is clamped around cylindrical jaws to reduce the gauge length to 3 cm. The tensile machine is piloted with an exponential head displacement up to L max = 15 cm to fix a constant strain rate ε . Photographs are shot during experiments to measure the deformation, and to detect a possible heterogeneity of deformation.

The experiments are done on the polystyrene at four temperatures (110°C, 120°C, 130°C, 140°C) and four strain rates for each temperature in the range 10 4 s 1 ≤ ε ≤ 10 2 s 1 . The temperature range was chosen to provide measurements in the rubber state above T g ≈ 90°C.

The deformation is not really homogeneous and an analysis of the local deformation was applied to measure the true stress as a function of true strain (Fig. 4). The rheology markedely changes with temperature. The stress and the strain hardening effect under large deformation strongly decrease with temperature. The strain rate sensitivity of the stress increases with temperature. This complex rheology can be represented by a one dimensional constitutive equation initially proposed by G'Sell and Jonas [9]. It has been selected in the context of polymer thermoforming:

σ = K p (T) 1 exp( wε) exp hε 2 ε m (1)
where T is the temperature, ε the von Mises equivalent strain rate, ε the von Mises equivalent strain, and K p , w, h and m are, respectively, the consistency, the viscoelastic strain hardening and strain rate coefficients of the polymer. These results are completed by data from literature [10,11] (Fig. 5). The consistency K p of the material decreases suddenly above T g and keeps low values (≈2 MPa.s m ) in the experimental domain. The so called viscoelastic coefficient w is found to be low (≈3.7) and independent on temperature. The low temperature behavior is characteristic of a solid like material, i.e., with a significant strain hardening and a low strain rate exponent (Fig. 5). At high temperature the polymer is a fluid like material, i.e., no strain hardening and a large strain rate coefficient. Consequently, the rheology changes in the temperature range 110°C 140°C, well above the glass transition. In this section, we will present briefly the mechanical approach used in the finite element model Tform 3, in which the previous constitutive equation has been implemented. The reader can refer to [12 13] for a more detailed description of the mechanical model and its numerical discretization. According to the membrane assumption, the deformed sheet is considered as a geometric surface, neglecting flexure and transverse shear. A material point is identified by two curvilinear coordinates: θ 1 and θ 2 which are, in the case of an initial sheet in the xy plane, of the initial x and y coordinate respectively. Hence, at any point M on the deformed sheet (vector ) the local tangent basis is defined by:

1 = Ž Žθ 1 2 = Ž Žθ 2 3 = 1 || 1 x 2 || 1 x 2 (2)
where x denotes the vectorial product. The equilibrium of the deformed sheet is expressed by the principle of virtual work (without inertia effects):

∀ ∀ ∀ ∀ * σ σ σ σ ij v i|j * ed Ω P 3 . * d Ω p . * d Ω c
(3) in which Einstein's convention of summation of repeated indices is used, as in the sequel. σ ij are the covariant components of the Cauchy stress tensor (plane stress: i,j = 1,2), . |j denotes covariant derivation with respect to the θ j material coordinate. In the first term, the virtual power dissipated by the normal and shear stresses are neglected compared to the virtual power dissipated by the normal and shear stresses is neglected compared the virtual power dissipated by in plane stresses due to the sheet expansion. Then the summation is extended only to i,j=1,2. Practically, a plane stress assumption is used, with only three components in the symmetric stress tensor. The variable e is the sheet thickness, are the friction stresses on regions Ω c on the sheet contacting the mold, P is the inflation differential pressure applied to the domain Ω p of the deformed surface Ω.

Starting with a balanced configuration Ω at time t, the problem consists of determining the unknown equilibrated configuration Ω' at t+∆t. Variables at t+∆t are denoted by "prime".

For application of eq. ( 3) at t+∆t, we have, for any velocity field :

σ σ σ σ' ij i|j * e'd Ω' P' ' 3 . * d Ω' p '. * d Ω' c (4)
This equation is solved for the incremental displacement field between Ω and Ω', provided that , e' and σ σ σ σ can be calculated from . Those relations are exposed hereunder with the resolution of (4).

• Contact and friction : the contact condition applied to Ω c may be either sticking (no relative velocity with respect to the mold) or sliding. In this case, the tangential stress is given by Coulomb's friction law (coefficient µ), in which the normal stress is the inflation pressure:

= µ P' (1 / | | | | ) (5)
• Thickness updating: the new local thickness e' is deduced from material incompressibility.

Denoting the metric tensor (g ij = . ), we have the following relation ( 6), which permits, for a given displacement field to evaluate and so update the thickness e' :

e' = e det ( ) / det ( )

• Resolution of incremental constitutive equations : the one dimensional constitutive equation ( 1) can be written as a classical viscoplastic power law:

σ = k ε m (7)
in which the k is the so called viscoplastic consistency. Hence, the flow rule derives from a viscoplastic potential Q. Under the assumption of isotropy, it yields :

ε ε ε ε = ŽQ Žσ σ σ σ = ŽQ Žσ Žσ Žσ σ σ σ = ε Žσ Žσ σ σ σ (8)
In convective curvilinear coordinates, the equivalent von Mises strain rate ε is expressed as :

ε 2 = 2/3 ε ij ε ij (9)
in which ε ij and ε ij denote the covariant and contravariant components of the strain rate tensor ε ε ε ε , respectively.

The equivalent von Mises stress has a similar expression, with the covariant and contravariant components of the deviatoric stress tensor :

σ = 3 2 s ij s ij = σ σ σ σ + p s ij = σ ij + p g ij (10)
Taking into account the above plane stress assumption, the hydrostatic pressure p is given by : p

= -1 3 Trace (σ) = -1 3 σ ij g ij i, j = 1,2 (11)
and the equivalent von Mises stress can be written :

σ = σ ij A ijkl σ kl (12) with A ijkl = 3 2 g ik g jl 1 2 g ij g kl i,j,k,l = 1,2 (13)
Using now the vectorial notation σ = σ 11 σ 22 σ 12 T , we have :

σ 2 = σ T σ (14)
where is a symmetric 3x3 matrix defined by :

A 11 =A 1111 A 22 =A 2222 A 33 =A 1212 A 12 =A 1122 A 13 =A 1112 A 23 =A 2212 (15)
The constitutive equation can then be written:

ε ε ε ε = ε σ σ σ σ σ = 1 k ε (1 m) σ σ σ σ (16)
where ε ε ε ε = ε 11 ε 22 ε 12 T . A semi implicit time integration scheme is used over the increment.

The incremental strain tensor ∆ε, the covariant components of which depend on the displacement according to (17) are written as (18):

∆ε ij = 1/2 (u i|j + u j|i + u m|i u |j m ) (17) ∆ε ε ε ε = ∆t (1 η) ε ε ε ε + η ε ε ε ε' (18)
Equations (16 18) clearly permit determination of the new local stress tensor σ', knowing the displacement field . A fully implicit scheme (η=1) is used (Euler backward scheme).

Finally, substitution of (5 6) and ( 16

= [ ] ε ε ε ε = [ ] (19)
The discretized form of the equilibrium equation can be written as :

( ) = [ '] T σ' e' d Ω' - P' [ ] T g 3 ' dS Ω p ' - [ ] T t' dS Ω c ' = 0 (20)
Hence, at every time step a non linear system for the vector of the nodal displacements is solved by the Newton Raphson method. It should be noted that a consistent tangent matrix ∂ /∂ can be expressed.

Regarding contact treatment, the simulation code can handle several tools (punch or molds)

that are supposed to be rigid. Their surfaces are discretized in triangular elements, so that it is possible to evaluate the relative positions of the nodes of the finite element mesh of the sheet, at each time increment. In case of slight penetration of a node in the tool, the normal displacement of this node comes back to the tangential tool surface. In addition, to prevent such undesirable effects, the nodes that are found in a thin boundary layer at the neighborhood of the tool surface are considered to be in contact and their normal displacement to the tool is prescribed equal to zero.

Thermal model and coupling

Considering the thinness of polymer sheets, the short processing times, and the low thermal diffusivity of polymers, it can be assumed that heat transfer is essentially one dimensional across the thickness of the sheet [14]. Consequently, s being the coordinate in the thickness direction (s=θ 3 ), the 1D heat transfer equation can be expressed:

ρc dT dt = Ž Žs λ ŽT Žs + σ σ σ σ:ε ε ε ε ( 21
)
where ρ is the specific mass, c the heat capacity and λ the heat conductivity. The following boundary conditions are accounted for. At sheet/air interface:

λ ŽT Žs sgn( ) = h conv (T T air ) ( 22
)
where h conv is the coefficient for heat exchange by convection, T air the air temperature and sgn( ) is ±1 depending on the orientation of the outward normal unit vector . At sheet/tool interface the surface temperature of the polymer sheet is prescribed by the interface temperature :

T inter = (b m T m +b s T s ) / (b m +b s ) ( 23
)
where b is the thermal effusivity λρc and T s is the average temperature in the thickness of the sheet. The initial temperature profile is assumed to be known at the beginning of the process. The following values of the thermal parameters characterize polystyrene: ρc = 1.92 10 6 J.m 3 .K 1 , λ = 0.2 W.m 1 K 1 , h = 33 W.m 2 .K 1 . Coulomb friction has been accounted for (µ = 0.4) [15,16].

Equation ( 14) is discretized in space and time, at each integration point of membrane, using a Galerkin 1D finite element method with linear three node elements (Fig. 6) and a semi implicit Crank Nicholson time integration scheme [14]. The coupling between the mechanical and the thermal solution is carried out at each time increment as explained in figure 7.

Multilayer formulation

The multilayer approach proposed here is based on the following assumptions:

• each layer is subjected to the same deformation as the mean surface of the sheet. This is consistent with the membrane approach and the thickness ratio of the different layers remains constant during the process.

• accordingly, the mechanical and thermal contacts are assumed perfect (no sliding, no thermal contact resistance) at interfaces between layers.

The algorithm is identical to the one for heat transfer coupling, except that all thermal and mechanical parameters used at integration points in thickness (ip_th) will now depend upon the material in which they are located (Fig. 8). Each material layer has an identified temperature dependent material law. A mean stress is then calculated from each stress contribution inside the thickness. Such a formulation is expected to be more precise than the reduction of the multilayer to a single "equivalent" material, especially when steep temperature gradients appear in the sheet thickness when one side of the multilayer contacts the tooling.

RESULTS

Axisymmetric thermoforming test

For the axisymmetric mold, only a 11°25' sector has been meshed with 160 triangles.

4.1.1 Thermal effect in the thickness of the sheet during the forming process

The initial orientation of the polystyrene sheet due to extrusion is relatively high, and it involves an evolution of the thickness during infrared heating as it is clamped. Consequently, the measured final thickness profile depends upon the measurement direction (extrusion or transverse direction, Fig. 9). The measurements shall be done along the extrusion direction to be consistent with the measured rheology.

The thermal effects shall be first considered. It is possible to predict the thermal gradient across the thickness of the sheet during the whole process (Fig. 10) and in the deeper part of the mold (Fig. 11). During the stretching in air the thermo mechanical model predicts a very low heating due to the deformation in the sheet and a very low cooling due to air convection at the surface. The deformation is nearly isothermal (temperature variation less than 10°C).

During this period the thickness of the sheet decreases (Fig. 10, 11). The polymer is quickly cooled after the contact with the mold and a large thermal gradient appears. A gradient of stress develops inside the thickness of the sheet as the result of the thermal dependence of the rheology. The contact on the insert is reached near the beginning of the forming process. The long heat transfer time inside the thick part allows a sliding of the sheet from 0.1s to 0.25s resulting in a continuation of the thinning of the film. After this period the cooling effect on the rheology is sufficient to balance in the plane stress and nearly halt the deformation (Fig. 10). On the contrary, the cooling of thinner film (higher deformation at the contact time) fixes almost instantaneously the polymer to the mold. The in plane stress is not sufficient to counterbalance thermal effects in the film. This mean the stress quickly increases in that position of the mean surface of the sheet. The deformation is nearly stopped in that position as an effect of the equilibrium relationship on the sheet (eq. 3). Consequently the thickness of the sheet in that position is nearly fixed by the cooling (Figs. 10 and11). This thermal effect coupled with a moderate but realistic friction coefficient (µ = 0.4) leads to a very efficient freezing of the polymer during contact with a low temperature mold (Fig. 12). The result is not so different from the sticking condition usually considered in most of the isothermal models. Sticking has no physical sense, since independent measurements always give a very low value for the friction coefficient of a solid polymer on a metal [15,16]. The present non isothermal formulation clearly uncouples the frictional and the thermal effects. If sticking is nevertheless considered there is no effect of the nature or temperature of the mold on the thickness profile and no possible processing action except the temperature control of the sheet. If the boundary condition is mainly dominated by the mold temperature this parameter can be used to favor a limited sliding to produce a more homogeneously thick part.

Prediction of processing time

The deformation of the sheet was simulated in various pressure histories and various temperatures for the polymer and the mold. The simulated time contact and measured values are systematically compared. When the differential pressure ∆P measured by the gauges is applied in the model the time of contact is always overestimated by a large factor, more than five. For example with the condition described in figure 10 the sheet is in contact at 0.095 s with the central insert when the model predicts more than 0.48 s to reach this position. The discrepancy between measurements and predictions is not explained. It must be noticed that the location of these pressure gauges give static pressure and should underestimate the pressure applied on the sheet.

Further calculations consider only the upper pressure P o . The predictions of time contact fit accurately the measurements for the various applied thermal and pressure conditions (Fig. 13). It is possible to analyze the effect of processing parameters on the time of contact either using the direct measurements or the predictions of the model. The upper pressure increases almost linearly with time up to the maximum pressure (Fig. 3). The measurements and the model show that the polymer is thermoformed during the transient period (Figs. 3 and 13).

The pressure rate shall be considered as the representative mechanical parameter and not the maximum pressure. It is the main mechanical parameter acting on the kinematics and hence on the time of contact (Fig. 14). An increase in the temperature of the polymer slightly decreases the time of contact (Fig. 14). This decrease is surprisingly low when the rheology shows a strong change of mechanical behavior in this temperature range (110 130°C) (Figs.

4 and 5). A higher polymer temperature leads to a higher strain rate coefficient (Fig. 5).

During the forming process the strain rate effect partly counterbalances the thermal effect on the consistency K p . These are the two main rheological parameters acting on the polymer deformation. The strain hardening coefficient h has a small effect in this range of deformation and is slightly more prominant near the end of the forming process. The mold temperature, up to T m = 60°C, does not significantly affects the time of contact, the slight change of the boundary condition induced by the cooling of the sheet in contact with the mold has a tiny effect on the remaining free part of the sheet.

Prediction of the final thickness

The thickness variation is calculated for each position throughout the thermoforming process. The local thickness can be shown along the radius of the part, more precisely on the projection along the radius. The mould geometry with a central insert is considered. The predicted and measured final thickness can be compared. The final thickness is calculated with various thermal (thermal or isothermal) and boundary (sticking or friction) conditions.

The thickness predicted by all the models is almost the same over most of the geometry. They differ significantly only at the contact with the insert. As it can be seen from figures 9 to 12, the isothermal model with friction (µ = 0.4) predicts an unrealistic thinning of the polymer on the insert (x ≤ 20 mm). This thinning comes from the sliding of the polymer sheet on the insert. The polymer rheology under isothermal condition is unable to halt the deformation. On the contrary a very large thickness is predicted on the insert when a sticking condition is applied with either an isothermal or non isothermal model. The thermal model is necessary to predict the thickness with areasonable value of the friction coefficient (µ = 0.4). It appears that thermal effects are very efficient at halting the polymer deformation. The deformation is then located in zones where the temperature is still higher, i.e., without any previous contact with the mold. In these thermal conditions, cooling has almost the same effect as the boundary conditions. The non isothermal model is much more representative of the local phenomena in thermoforming. It is able to account independently for interface tribology and heat transfer, and finally results in a very good agreement between measured and computed thickness profiles.

Orientation in polystyrene part

Model and experiments predict the same heterogeneous deformation. The thinning, mostly located in the corners, results from an elongation which can induce an orientation of the polymer. This orientation acts to increase the mechanical properties. It can partly counterbalance the local thinning of the industrial part. The size and geometry of the part lead to difficult measurements of the local mechanical properties. The local polymer orientation , measured by birefringence, was used to check the possible heterogeneity of the mechanical properties. A K compensator on a Reichert Zetopan pol microscope was used to measure the optical path difference ∆λ. The birefringence is deduced from the local thickness knowledge e, ∆n = ∆λ /e. These birefringence measurements were done along the radius for two processing conditions (Fig. 15). The strong condition (high pressure P = 0.4 Mpa, T s = 130°C) minimizes the relaxation of polymer chains and leads to a higher orientation than the soft condition (low pressure P = 0.1 Mpa, T s = 140°C). Both conditions produce the highest orientation in the inner corner of the insert (contact sensor C 4 on figure 2) and the lowest at the first contact with the insert (contact sensor C 1 ). The very low orientaion on the insert results from the very low deformation necessary to reach the contact.No significant sliding appears during the contact of the polymer on this part of the mould. Consequently, no processing effectcan be found in this part of the polymer. On the contrary, the highest deformation (contact sensor C 4 ) is the most sensitive to processing parameters. The strong conditions must induce the higher orientation and thus the higher mechanical properties.

Multilayer application

We have simulated the forming of a bilayered sheet : 80% polystyrene (PS), 20% polypropylene (PP), total thickness 1 mm. The rheological parameters for both polymers are those given in section 2. The sheet is formed successively with a plug with a 50 mm radius (maximum ram displacement ∆z= 30 mm at 0.15 s) and pressure (linear increase of 0.8 MPa between 0.15 s and 0.65 s). The initial temperatures are 20°C for tools and 150°C for the sheet.

The generic shape of the thickness distribution along the radius is due to the punch assisted forming. It is shown in table 1 and Figure 16 that the deformation of the bi layer is not a simple addition or average of the deformation of each layer component with the same thickness. In addition, the results clearly depend on which material is on top. This is due to the fact that the cooling effect is different for polystyrene or polypropylene : the consistency of polystyrene decreases suddenly above T g whereas it follows an Arrhenius law for polypropylene. Also m and h are temperature dependent for polystyrene but constant for polypropylene. This point is to be confirmed experimentally.

Non axisymmetrical model

In order to test the robustness and the results of the code, the forming of a one layer shallow component for food packaging has been studied (Fig. 17). Only a quarter of the symmetric industrial forming mold has been computed, using 7591 and 6463 triangles for the sheet and mold meshes respectively. The material is polypropylene, initially 0.475 mm thick.

The following parameters can be considered realistic for industrial forming conditions :

aluminum mold (5°C), linear pressure reaching 0.7 MPa at 0.5s, air temperature 20°C, initial sheet temperature 150°C. As shown on Fig. 18, the central bulging of the sheet is very fast, reaching the bottom at 0.04s. The forming is then slowed down until the end. Figure 19 shows the deformed finite element mesh at the end of the process. Iso values of temperature (not shown here) indicate that during forming the free regions are the warmer (close to the initial temperature), and that in contacting regions, the longer the contact duration, the cooler is the four values per point. The experimental dispersion has been found high (up to 18% of the average value). The average of the relative errors between computed and average experimental thickness is reasonable : 10% at the bottom, 18% at the wall.

CONCLUSION

A thermo mechanical model including thermal evolution during the process is necessary to predict the deformation, and the final thickness of thermoformed parts. The thermal effect during contact with the tools mainly controls the sliding of the material on the mold. The model includes a temperature dependent visco plastic rheology whose parameters are measured for a polystyrene under uniaxial elongation with constant strain rate and isothermal conditions. The polystyrene exhibits a rheological change from liquid to solid like behaviour over a wide temperature range.

ACKNOWLEDGMENT

The experimental part of this work has been done at Elf Atochem (Cerdato) which has also supported this study. In the Elf Atochem company, Daniel Lebouvier, Jean Claude Jammet and Ronan Andro must be thanked for their support, Stéphane Galliot for the instrumentation and his help for on line measurements. Prof. J.M. Haudin is thanked for helpful discussions on this paper. 1D finite element mesh membrane finite element mesh

Figure 6 : One dimensional approach through the sheet thickness for thermal coupling.

A/ Thermal resolution • for all membrane elements • for all integration points of the element perform a 1D finite element computation of the new temperature profile T' across the thickness. The source term is given by the mechanical resolution of previous time step.

for all integration points "ip_th" of the 1D mesh update the temperature dependent coefficients: K' p (T'), m'(T') and h'(T').

B/ Mechanical resolution B1) computation of the residual vector ( )

• for all membrane elements • for all integration points of the element for all integration points "ip_th" of the 1D mesh solve the constitutive equations for σ .27 .23 .22 .22 .22 .21 .21 .21 .20 .19 .19 .18 .20 .20 .17 .19 .21

.

23

.

26

.

29

.33

.36

.30

.

42

.47

.CALCULATED THICKNESS (mm)

EXPERIMENTAL THICKNESS (mm)

.

.214 .258 .252

-40.

-30.

-20.

-10. "diag2" directions on figure 19).

  This rheological model is able to describe the mechanical behavior of a great number of polymers, either amorphous or semi crystalline, for a large temperature range. The parameters for polypropylene have been identified for a large temperature range (20°C ≤ T ≤ 150°C) by means of tensile tests on sheets with the same range of deformation rate [9] : K p = 5.31 10 3 exp [2.85 10 3 / T] MPa.s m w = 30, h = 0.4, m = 0.04 For polypropylene only the consistency K p significantly depends on temperature, the value of the other parameters remain almost constant. These rheologies, model and parameters, are included in the finite element code to predict the polymer deformation coupled with thermal effects. 3 3 D THERMO MECHANICAL MODEL FOR A MULTILAYERED MATERIAL 3.1 Mechanical model and resolution

  18) in the equilibrium equation (4) leads to a non linear equation for the displacement field . Its spatial discretization by the standard Galerkin finite element method (linear triangles or quadrangles) is detailed by Bellet [12, 13]. Let us denote [ ] the interpolation matrix and [ ] the differentiation matrix defined by :

  Figure 20. Experimental points are issued from measurements on two different parts, yielding

Figure 4 :Figure 5 :

 45 Figure4: Uniaxial stress strain curves measured under constant strain rate (5 10 3 s 1 ) in isothermal conditions for four temperatures.

Figure 7 :Figure 8 :

 78 Figure 7 : Thermo mechanical coupling algorithm.

Figure 9 :Figure 12 :

 912 Figure 9 : Comparison between computed (solid line) and measured thickness distribution along the extrusion and transverse directions after an axisymmetrical forming with a central insert.

Figure 20 :

 20 Figure20: Example of comparison of measured and computed thickness a: "length" and b: "diag2" directions on figure19).

Table 1 -

 1 Computed thickness profile of bi layered films (polypropylene or polystyrene on top) formed with a punch compared with the pure polymers: Temperature = 150 °C; Pressure rate = 1.6 MPa/s. Prediction of forming time for the bilayer and pure films (see processing conditions in the text).Figure 19 : Final deformed finite element mesh and selected directions for the thickness comparison.

	R(mm)


Marie Hélène Vantal

Bernard Monasse

Michel Bellet

Numerical simulation of the thermoforming of multi layer polymer sheets

This paper presents a 3D finite element model of the thermoforming process. The polymer material is assumed to obey the viscoplastic law proposed by G'Sell and Jonas, the coefficient of which have been identified for a polystyrene. The implementation of this constitutive model, the heat transfer coupling and the proposed multi layer approach are detailed in the paper. Axisymmetric validation tests and an application to the thermoforming of an industrial component are also reported.

INTRODUCTION

The thermoforming process consists in heating a polymer sheet and shaping it inside a mold. The deformation results from a pressure cycle eventually coupled with the use of a moving tool. The forming of thin products made of various polymer materials is carried out by this method: rubbery amorphous, solid semi crystalline or multi layer composites. The heating temperature depends upon the selected polymer. At low temperature, the forming is limited by the too high rigidity of the sheet. On the contrary, at high temperature, the sheet deforms by gravity and the forming is very difficult to control. Hence, the thermoforming is operated above the glass transition temperature (T g ), in the rubber like behavior domain for amorphous polymers, and close to the melting temperature for semi crystalline polymers.

The main problem of the process is the thinning in the corners of the parts, which leads to a decrease of the mechanical properties of the shaped components. The optimization of the final thickness profile is generally achieved by trial and error, changing the design of the component, the polymer material, and the process parameters such as the heating temperature distribution, the mold temperature, the pressure curve, or using a punch for deep parts. The numerical modelling should then result in a more efficient optimization of the process.

Various numerical simulation models have been proposed [START_REF] Warby | Finite element model of viscoelastic membrane deformation[END_REF][START_REF] Lorenzi | Finite element simulation of thermoforming and blow molding[END_REF][START_REF] Kouba | Computer simulation of thermoforming in complex shapes[END_REF][START_REF] Shrivastava | Large deformation finite element analysis of viscoelastic membranes with reference to thermoforming[END_REF]. Regarding the application to three dimensional formings, they are generally based on the membrane mechanical approximation, associated with the finite element method. Those computations use either hyperelastic [START_REF] Lorenzi | Finite element simulation of thermoforming and blow molding[END_REF][START_REF] Shrivastava | Large deformation finite element analysis of viscoelastic membranes with reference to thermoforming[END_REF] or viscoelastic [START_REF] Warby | Finite element model of viscoelastic membrane deformation[END_REF][START_REF] Kouba | Computer simulation of thermoforming in complex shapes[END_REF]) constitutive equations, without any heat transfer coupling and a sticking contact hypothesis at polymer tools interface.

It should be pointed out that this latter assumption has no experimental evidence. Most of the time, it is used just because the computation is isothermal and cannot account for the decrease of the polymer temperature and its "freezing" after tool contact. Another thermomechanical coupling is the high self heating source term due to the high strain rates. As the behavior of polymers is known to be highly temperature dependent, it seems essential to couple heat transfer and mechanical models.

In the present paper, we describe a 3D finite element model with membrane approximation to predict the deformation. Concerning the contact conditions with punch and mold, either sliding contact with Coulomb's friction law or sticking contact can be considered. The paper is focused on the thermal evolution effects during the forming, and on the approach to multilayer sheet forming.

MATERIAL BEHAVIOR IDENTIFICATION

The one dimensional constitutive equation initially proposed by G'Sell and Jonas (1979) has been selected in the context of polymer thermoforming: σ = K p (T) 1 exp( wε) exp hε 2 ε m (1) where T is the temperature, ε the von Mises equivalent strain rate, ε the von Mises equivalent strain. K p , w, h and m are material parameters.

This law accounts for the material behavior of a great number of polymers, either amorphous or semi crystalline, in a large temperature interval. [START_REF] Duffo | Rheology of polypropylene in the solid state[END_REF] have identified parameters for polypropylene by means of tensile tests on sheets, yielding: K p = 5.31 10 3 exp [2.85 10 3 / T] MPa.s m w = 20, h = 0.4, m = 0.087

We have studied the material behavior of a polystyrene compound, by uniaxial tensile testing at constant axial strain rate. Tests have been carried out at four constant temperatures: 110, 120, 130 and 140°C (which cover the forming temperature domain of this polymer) and five different strain rates between 10 4 and 10 2 s 1 . The associated determination of the coefficients of equation (1) shows the transition from a solid like behavior to a liquid like one: see figure 1. The consistency K p of the material decreases suddenly above T g and keeps low values (≈2) in the experimental domain. The so called viscoelastic coefficient w is found low (≈3.7) and independent on temperature. Finally, it has appeared that the basic equation (1) has to be modified, in order to include a temperature dependency of the strain rate sensitivity coefficient m and the hardening coefficient h.

MECHANICAL FORMULATION

" # $ #
According to membrane mechanical assumption, the deformed sheet is considered as a geometric surface, neglecting flexion and transverse shear. A material point is identified by two curvilinear coordinates: θ 1 and θ 2 which are, in the case of an initial sheet in xy plane, the initial x and y coordinate respectively. At any point of the deformed sheet (vector ) the local tangent basis is defined by:

1 = Ž Žθ 1 2 = Ž Žθ 2 3 = 1 || 1 ∧ 2 || 1 ∧ 2 (2)
The equilibrium of the deformed sheet is expressed by the principle of virtual work (without inertia):

∀ * σ σ σ σ ij v i|j * ed Ω P 3 . * d Ω p . * d Ω c = 0 (3)
where σ ij are the covariant components of the Cauchy stress tensor (plane stress: i,j=1,2), |j denotes covariant derivation with respect to θ j , e is the sheet thickness, are the friction stresses on regions Ω c contacting the mold, P is the inflation pressure applied to the domain Ω p of the deformed surface Ω.

"% & ' " (

Starting with a balanced configuration Ω at time t, the problem consists in determining the unknown equilibrated configuration Ω' at t+∆t. Variables at t+∆t are denoted "prime". By application of eq.( 3) at t+∆t, we have, for any velocity field *:

σ σ σ σ' ij i|j * e'd Ω ' P' ' 3 . * d Ω 'p '. * d Ω 'c = 0 (4)
This equation is solved for the incremental displacement field between Ω and Ω', provided that ', e' and σ σ σ σ' can be calculated from . Those relations are exposed hereunder and the resolution in 3.3.

• ) * as regard the contact condition applied to Ω c , it may be either sticking (no relative velocity with respect to the mold) or sliding: in this case, the tangential stress is supposed to be given by the Coulomb's friction law ( coefficient µ), in which the normal stress is the inflation pressure:

' = µ P' (1 / || ||) (5)
• &! + ,* the new local thickness e' is deduced from material incompressibility. Denoting the metric tensor (g ij = i . j ), we have:

e' = e det( )/det( ') (6) • ) - $
* the one dimensional constitutive equation (1) can be written as a classical viscoplastic power law:

σ = k ε m (7)
Hence, the flow rule derives from a viscoplastic potential Q. Under isotropy assumption, it yields:

ε ε ε ε = ŽQ Žσ σ σ σ = ŽQ Žσ Žσ Žσ σ σ σ = ε Žσ Žσ σ σ σ (8)
In convective curvilinear coordinates, we have:

ε 2 = 2/3 ε ij ε ij σ 2 = 3/2 s ij s ij = σ σ σ σ T σ σ σ σ (9) A ijkl = 3/2 g ik g jl 1/2 g ij g kl (10)
The constitutive equation can then be written:

ε ε ε ε = (ε/σ) σ σ σ σ = (1/k) ε (1 m)
σ σ σ σ (11) A semi implicit time integration scheme is used over the increment. The incremental strain tensor ∆ε ε ε ε, the covariant components of which depend on the displacement according to ( 12) is written as (13):

∆ε ij = 1/2 (u i|j + u j|i + u m|i u |j m ) (12) ∆ε ε ε ε = ∆t (1 η) ε ε ε ε + η ε ε ε ε' (13
) Equations (11 13) clearly permit to deduce the new local stress tensor σ σ σ σ', knowing the displacement field . Practically, a fully implicit scheme (η=1) is used.

" * ) ! .)" " "/
Finally, injection of (5 6) and (11 13) in the equilibrium equation ( 4) leads to a non linear equation for the displacement field . Its spatial discretization by f.e.m. (linear triangles or quadrangles) is detailed by [START_REF] Bellet | Modélisation numérique du formage superplastique de tôles[END_REF][START_REF] Bellet | Numerical simulation of thin sheet forming processes by the finite element method[END_REF]. At every time step a non linear system for the vector of the nodal displacements is solved by the Newton Raphson method, with a consistent tangent matrix.

HEAT TRANSFER RESOLUTION

Considering the thinness of polymer sheets, the short processing times, and the low diffusivity of polymers, it can be assumed that heat transfer is essentially one dimensional across the thickness of the sheet [START_REF] Vantal | Etude numérique et expérimentale du thermoformage des polymères[END_REF]). Consequently, s being the coordinate in the thickness direction (s=θ 3 ), the 1D heat transfer equation can be expressed:

ρc dT dt = Ž Žs λ ŽT Žs + σ σ σ σ:ε ε ε ε (14)
where ρ is the specific mass, c the heat capacity and λ the heat conductivity. The following boundary conditions are accounted for. At sheet/air interface: λ ŽT Žs sgn( ) = h conv (T T air ) (15) where h conv is the coefficient for heat exchange by convection, T air the air temperature and sgn( ) is ±1 depending on the orientation of the outward normal unit vector . At sheet/tool interface, due to the much higher diffusivity of metals, we will assume that the surface temperature of the polymer sheet is prescribed to the interface temperature given by: T inter = (b tool T tool +b sheet T sheet ) / (b tool +b sheet ) (16) where b is the thermal effusivity λρc and T sheet is the average temperature of the sheet in the thickness. The initial temperature profile is assumed to be known at the beginning of the process.

Equation ( 14) is discretized in space and time, at each integration point of membrane finite elements, using a Galerkin 1D f.e.m. (figure 2) and a semi implicit time integration scheme [START_REF] Vantal | Etude numérique et expérimentale du thermoformage des polymères[END_REF].

The coupling between the mechanical and the thermal resolution is carried out at each time increment as explained in figure 3.

MULTI LAYER FORMULATION

The multi layer approach proposed here is based on the following assumptions:

• each layer is submitted to the same deformation as the mean surface of the sheet. This is consistent with 1D finite element mesh membrane finite element mesh Fig. 2 1D approach for thermal coupling A/ Thermal resolution • for all membrane elements • for all integration points of the element perform a 1D f.e. computation of the new temperature profile T' across the thickness. The source term is given by the mechanical resolution of previous time step.

for all integration points ipth of 1D mesh update the temperature dependent coefficients: K' p (T'), m'(T') and h'(T').

B/ Mechanical resolution B1) computation of the residual vector ( ) • for all membrane elements • for all integration points of the element for all integration points ipth of 1D mesh solve the constitutive equations for σ the membrane approach in which transverse shear and flexion are neglected. Hence, the thickness ratio of the different layers remains constant during the process.

• accordingly, the mechanical and thermal contact are assumed perfect (no sliding, no thermal contact resistance) at interfaces between layers.

The algorithm is identical to the one for heat transfer coupling, except that all thermal and mechanical parameters used at integration points in thickness (ipth) will now depend upon the material in which they are located (see fig. 4). Each layer material has an identified temperature dependent behavior law.

Such a formulation is expected to be more precise than the reduction of the multi layer to a single "equivalent" material, especially when steep temperature gradients appear in the sheet thickness when one side of the multi layer contacts the tooling. The mold geometry is axisymmetric (diameter ≈ 140 mm, depth 60 mm), including a central insert in the bottom (see fig. 6) which makes it very sensitive both to cooling and friction effects. The forming parameters are the followings: aluminium mold, 20°C, air temperature = 20°C, linear pressure vs time: 0.6 MPa.s 1 , initial sheet temperature 130°C. The one layer polystyrene sheet is initially 1.1 mm thick. Its orientation due to extrusion is relatively high, and it involves an evolution of the thickness during infrared heating as it is clamped. Consequently, the measured final thickness profile depends upon the measurement direction (extrusion or transverse direction, see fig. 7).

0"% ) !

Due to axisymmetry, only a narrow sector has been meshed with 160 triangles. We have used the rheological data of fig. 1 and the following values of thermal parameters: ρc = 1.92 10 6 J.m 3 .K 1 , λ = .2 W.m 1 .K 1 , h = 33 W.m 2 .K 1 . A Coulomb friction has been accounted for (µ=0.4).

The predicted final thickness profile, using the coupled thermomechanical model, is in good agreement with the measurements (figure 7).

Regarding kinematics, it should be noticed that the measured differential pressure was found very low. The use of these values as a prescribed boundary condition (P(t) in equation ( 3)) has lead to much slower forming rates than the actual ones. However, if the measured upper pressure is used instead of the differential pressure, then computed forming times are in agreement with experimental ones, as shown in table 1. This unexpected result needs more investigation. contact time (s) at sensor # 1 2 3 experiment 0.095 s 0.400 s 0.510 s computation 0.102 s 0.311 s 0.390 s Table 1 0" ( ) ) ! ,

On figure 8, three different computed thickness profiles have been plotted, using three different computational options.

As already said, the assumption of sticking contact has no clear experimental evidence. However, as shown on figure 8, it is the only means to get reasonable results if an isothermal model is used. For example the combination (isothermal ; high friction) yields completely erroneous thicknesses (far too low on the insert, despite the high friction coefficient)! On another hand, the use of the present non isothermal formulation permits to decouple clearly the frictional and the thermal effects: the quick cooling of the polymer when contacting the mold increases the consistency and the hardening coefficient of the material for the contact zones, as the. strain rate sensitivity coefficient decreases. The evolution of these coefficients localize the deformation in the warmer zones: the strain rate values are almost zero in the "frozen" zones. The non isothermal model is much more representative of the local phenomena in thermoforming and is able to account independently for interface tribology and heat transfer. Fig. 8 Computed thickness for three models.

0" 2 1

We have simulated the forming of a bi layer: 80% polystyrene (PS), 20% polypropylene (PP), total thickness 1mm. The rheological parameters for both polymers are those given in section 2. The sheet is formed successively with a punch (z= 30mm at 0.15 s) and pressure (linear increase of 0.8 MPa between 0.15 s and 0.65 s). The initial temperatures are 20°C for tools and 150°C for the sheet.

It is shown in with the same thickness. In addition, the results clearly depend on which material is on top. This is due to the fact that the cooling effect is different for PS or PP: the consistency of PS decreases suddenly above T g whereas it follows an Arrhenius law for PP. Also m and h are temperature dependent for PS whereas constant for PP.

APPLICATION: INDUSTRIAL FORMING

In order to test the robustness and the results of the code , the forming of a one layer shallow component for food packaging has been studied (fig. 10). The material is polypropylene, initially 0.475 mm thick. The industrial forming conditions have not been accurately measured, but the following figures can be considered realistic: aluminum mold (5°C), linear pressure reaching 0.7MPa at 0.5s, air temperature 20°C, initial sheet temperature 150°C.

Only a quarter of the symmetrized actual forming has been computed, using 7591 and 6463 triangles for the sheet and mold meshes respectively. A Coulomb's friction law (^=0.4) is assumed at the mold surface, except near the edge where sticking contact is prescribed. As shown on fig. 12, the central bulging of the sheet is very fast, reaching the bottom at 0.04s. The forming is then slowed down until the end. Figure 11 shows the deformed finite element mesh at the end of the process. Iso values of temperature (not shown here) indicate that during forming the free regions are the warmer (close to the initial temperature), and that in contacting regions, the longer the contact duration, the cooler is the sheet.

The final thickness of actual parts has been measured along the five directions mentioned on figure 11 -30.

-20.

-10. 0.0 10.

.CALCULATED THICKNESS (mm) EXPERIMENTAL THICKNESS (mm) Fig. 13 Example of comparison between measured and computed thickness ("length" on figure 12). average value). The average of relative errors between computed and average experimental thickness is good: 10% at the bottom, 18% at wall.

CONCLUSION

Thermomechanical coupling and a multi layer approach have been implemented in a finite element membrane model. They are shown to be very efficient tools to improve the predictive character of finite element simulations of thermoforming.

ACKNOWLEDGEMENT

The experimental part of this work has been done at Elf Atochem (Cerdato) which has also supported this study.

Fig. 1

 1 Fig.1 Coefficients m, h and K p (MPa.s m ) of equ. (1) for polystyrene vs temperature. (Grey dots identified by tensile tests, black squares issued from literature).

  Fig.3 Thermo mechanical coupling algorithm.

Fig. 6

 6 Fig.4 Multi layer formulation (schematic).6 TEST OF THE MODEL

  Fig. 7 Comparison between measured and computed thickness distribution.

Fig. 10

 10 Fig.10 Selected polypropylene test part.

  Fig.11 Final deformed finite element mesh and selected directions for thickness comparison.

Fig. 12

 12 Fig. 12 Computed deformed sheet

table 2

 2

			and figure 9 that the
	deformation of the bilayer	a simple addition or
	average of the deformation of each layer component
	PP only	PP_PS	PS_PP	PS only
		PP on top	PS on top
	0.59 s	0.50 s	0.39 s	0.45 s

Table 2

 2 Predicted forming times.

	R(mm)

  . An example of comparison with the computed values is shown in figure 13. Experimental points are issued from measurements on two different parts, yielding four values per point. The experimental dispersion has been found high (up to 18% of the

	Z(mm)
									.48
									.30
								.398	.25
									.26
								.36	.25
									.25
							.25	.23
	.325	.34			.30	.27 .23	.18	.22
	.35	.34	.32 .29 .27	.26 .27	.30	.30	.28 .27	.22 .21
	.0 -40.		18.6	37.2		55.8	length (mm) 74.4 93.0


Michel Mollard
email: michel.mollard@ujf-grenoble.fr

On Vertex Partitions of Hypercubes by Isometric Trees

Keywords: Graph, Perfect code, Hypercube, Vertex Partition, Tiling

When n = 2 m -1 M.Ramras proved, by a counting argument, that for any isometrically embedded tree T on n edges in Q n there exists a group of translations G such that {g(T ); g ∈ G} is a vertex partition of Q n . Considering a more general context we are able to give an explicit construction of G and can construct non group vertex partitions by isometric trees. We extend also this problem to vertex partition of Q n by translates of an isometrically embedded tree on n = 2 m -1 edges for any n ≥ n.

Introduction

Twenty years ago M.Ramras [START_REF] Ramras | Symmetric Vertex Partitions of Hypercubes by Isometric Trees[END_REF] published a paper where he answered the following question of D. Rogers: If n = 2 m -1 does the hypercube Q n have a vertex partition into antipodal paths? M.Ramras gave an explicit construction of such a partition. For this purpose he exhibited a set of generators of a subgroup G of the group of translations Σ(Q n ) ⊂ Aut(Q n ), such that the set of translates {g(P ); g ∈ G} of the vertex set of P of an antipodal path is the desired vertex partition.

He extended this result, proving, via a counting argument, that for any tree on n = 2 m -1 edges isometrically embedded in Q n with vertex set P there exists a subgroup G of Σ(Q n ) such that the set of translates{g(P ); g ∈ G} is a vertex partition. Notice that, in the general case, the author's method does not give an explicit construction of G.

This nice result generalizes the existence of perfect single-error-correcting codes constructed first by R.W.Hamming [START_REF] Hamming | Error detecting and error correcting codes[END_REF]. In this case, we take as tree the star K 1,n . Notice that J.L.Vasiliev [START_REF] Vasilev | On ungrouped, close-packed codes (in Russian)[END_REF] constructed, for n = 2 m -1, n ≥ 15, perfect codes which are not equivalent to linear codes, i.e. vertex partitions by stars such that the set of translations G is not a subgroup, or the translate of a subgroup, of the group of translations Σ(Q n ).

In another series of papers M.Ramras considers edges partitions of Q n into isomorphic trees. More recent work have been done of this subject but is seems that this is not the case for vertex partitions. Both problems arise in the context of parallel computing and thus it will be interesting to improve our knowledge of vertex partitions.

Our goal at the beginning of this work was to prove the existence of non group partitions into antipodal paths. It seems also interesting, for the general case of trees, to give an explicit construction of a group, or more generally a set, of translations G. It would be also nice, for the case of the path P n to understand how the group proposed by M.Ramras can be derived from the Hamming code. In fact we found that all these problems can be easily solved, using elementary linear algebra.

We arrive at the conclusion that looking for a vertex partition of Q n by translates of an isometrically embedded tree on n edges is a problem independent, in some sense, of the choice of the tree, thus is equivalent to looking for a perfect code.

Definitions and main result

Let F n be the vector space of dimension n over the finite field Z 2 . The hypercube of dimension n is the graph Q n whose vertices are the vectors of F n , and where two vertices are adjacent if they differ in exactly one coordinate.

The Hamming distance between two vectors x, y ∈ F n , d(x, y) is the number of coordinates in which they differ. Notice that Hamming distance is the usual graph distance on Q n .

The support of a vector x is the set {i ∈ {1, 2, . . . n} ; x i = 0}. The parity function is the function from F n to Z 2 defined by π(x 1 , x 2 , . . . , x n ) = x 1 + x 2 + . . . + x n .

A perfect code, or more precisely a perfect single-error-correcting code of length

n is a set C of vertices of Q n such that every vertex x ∈ V (Q n ) is at distance at most 1 of exactly one element of C.
Two codes C and C are called equivalent if C can be obtained from C by an automorphism of Q n , thus by applying a translation from a fixed vector and a fixed permutation of the coordinates. Using a translation of all the perfect code vectors by one of them we can always assume that the zero vector 0 belongs to the code.

Let e 1 , e 2 , . . . , e n be the standard basis of F n , thus e i denotes the vector with just one single non zero coordinate position i. Let e 0 = 0 be the zero vector and let 1 = (11 . . . 1). Denote by B n the set {e 0 , e 1 , e 2 , . . . , e n }. The direction of an edge

xy of Q n is the integer i ∈ {1, 2, . . . , n, } such that y = x + e i .
For a subset A of vectors of F n and a vertex x, let x+A be the set {x + a; a ∈ A}.

By definition of a perfect code

C = {c 1 , c 2 , . . . , c k } is a perfect code if and only if the sets c 1 + B n ,c 2 + B n ,. . . ,c k + B n define a partition of F n i.e. • F n = c 1 + B n ∪ c 2 + B n ∪ . . . ∪ c k + B n • ∀i, j ∈ {1, 2, . . . , k} c i + B n ∩ c j + B n = ∅
This immediately leads to a necessary condition for the existence of C, the so-called packing condition, 2 n = |C|(n + 1) thus n = 2 m -1 for some m.

By analogy, G.Cohen, S.Litsyn, A.Vardy and G.Zémor [START_REF] Cohen | Tilings of Binary spaces[END_REF] define a set S as a tile of F n if there exists a set C = {c 1 , c 2 , . . . , c k } such that the sets c 1 +S,c 2 +S,. . . ,c k +S define a partition of F n . Notice that the definition is symmetric, C is also a tile of F n , and they call the pair(C, S) a tiling of F n . R.W.Hamming [START_REF] Hamming | Error detecting and error correcting codes[END_REF] constructed, for any integer m, a linear subspace of F n , where n = 2 m -1, which is a perfect code . It is easy to prove that all linear perfect codes are Hamming codes. In 1961 J.L.Vasiliev [START_REF] Vasilev | On ungrouped, close-packed codes (in Russian)[END_REF], and later many authors ( [START_REF] Cohen | Covering Codes, Chap. 11[END_REF][START_REF] Heden | A survey on perfect codes[END_REF] for a survey) constructed perfect codes which are not linear codes.

Let W = (v 1 , v 2 , . . . , v n ) be a basis of F n . We will denote by θ W be the auto-

morphism of F n defined by θ W ( n i=1 λ i e i ) = ( n i=1 λ i v i ). Lemma 1 Let V = {v 0 , v 1 , . . . , v n } where W = (v 1 , v 2 , . . . , v n ) is a a basis of F n and v 0 = 0. Then C = {c 1 , c 2 , . . . , c k } is a perfect code if and only if θ W (c 1 ) + V, θ W (c 2 ) + V, . . . , θ W (c k ) + V is a partition of F n .
Proof : Notice that for any j ∈ {0, 1, . . . , n} we have θ W (e j ) = v j . By linearity for any i, j we have θ

W (c i ) + v j = θ W (c i ) + θ W (e j ) = θ W (c i + e j ). Therefore, because θ W is an automorphism, θ W (c i ) + v j = θ W (c i ) + v j if and only if c i + e j = c i + e j .
Assume that C is a perfect code then the θ W (c i ) + V are disjoint. Furthermore for any x of F n we knows that there exist i ∈ {1, 2, . . . , k} and j ∈ {0, 1, . . . , n} such that θ -1 W (x) = c i + e j , thus x = θ W (c i ) + v j and we have a partition.

Conversely if θ W (c 1 ) + V, θ W (c 2 ) + V, . . . , θ W (c k ) + V is a partition of F n then the θ W (c i ) + B n are disjoint.
Moreover for any x of F n there exist i and j such that

θ W (x) = θ W (c i ) + v j and thus x = c i + e j . 2
For two graphs G and H an isometric embedding of G in H is a map α : V (G) →

V (H) which preserves distance. By extension we will denote by α(G) the subgraph of H induced by α(V (G)). If G is injectively embedded in Q n we will say that there exists a vertex partition of Q n by G if there exists a tiling of V (Q n ) by α(V (G)).

It is immediate to check, as noticed by M.Ramras [START_REF] Ramras | Symmetric Vertex Partitions of Hypercubes by Isometric Trees[END_REF], that a tree T is isometrically embedded in Q n if and only if no edges of α(T ) use the same direction. If α is an isometric embedding in Q n then for any translation t the map α = α + t is also an isometric embedding. Therefore, if a graph G is isometrically embeddable in Q n then, for any vertex x of G there exists an isometric embedding such that α(x) = 0.

Lemma 2 Let T be any tree on p ≤ n edges, and let α be an isometric embedding

of T in Q n . Assume α(T ) = {0, v 1 , v 2 , . . . , v p } then the vectors v 1 , v 2 , . . . , v p are linearly independent.
Proof : The proof is by induction on p. The result is clearly true when p = 1 and assume it holds for any tree on p -1 edges. Let x be a terminal vertex of T and let xy be the edge of T incident to x. Consider the tree T obtained by deletion of x from T . We can always assume that α(x) = 0 thus α(x) = v i for some i ∈ {1, 2, . . . , p}. Let j ∈ {1, 2, . . . , n} such that α(x) = α(y) + e j . The restriction to T of α is an isometric embedding thus the vectors {v k ; k ∈ {1, 2, . . . , p} , k = i} are independent by induction hypothesis. Notice that e j does not not appear in the basis decomposition of the {v k ; k ∈ {1, 2, . . . , p} , k = i}. But α(y) belongs to this set and because v i = α(y) + e j the vector v i is also linearly independent of them.

2

Theorem 3 Let T be any tree on n = 2 m -1 edges, and let α be an isometric embed- Proof : By lemma 2 W is a basis of F n and the result follows by lemma 1.

ding of T in Q n . Assume α(T ) = {0, v 1 , v 2 , . . . , v n }. Then the vectors v 1 , v 2 , . . . , v n form a basis W of F n . Furthermore if C = {c 1 , c 2 , . . . , c k } is a perfect code of Q n then {θ W (c 1 ), θ W (c 2 ), . . . , θ W (c k )}
Conversely if S = {s 1 , s 2 , . . . , s k } defines a vertex partition of Q n by T then

θ -1 W (s 1 ), θ -1 W (s 2 ), . . . , θ -1 W (s k
) is a perfect code and thus all vertex partition of Q n by T arise in this way.

2

Notice that the set

D = {θ W (c 1 ), θ W (c 2 ), . . . , θ W (c k )} is a linear subspace if and only if C is linear. Furthermore if b 1 , b 2 . . . , b p is a basis of C then θ W (b 1 ), θ W (b 2 ), . . . , θ W (b p )
will be a basis of D.

Consider now a vertex partition of Q n by translates of an isometrically embedded tree on n edges for some n ≥ n. By the packing condition, 2 n = (n + 1)|D|, thus n = 2 m -1 for some m.

Corollary 4 Let T be any tree on n = 2 m -1 edges, and let α be an isometric

embedding of T in Q n , n ≥ n.
Then there exits a vertex partition of Q n by translates of the embedded tree α(T ).

Proof : The vertices of α(T ) = {0, v 1 , v 2 , . . . , v n } define a subspace of dimension n. By a permutation of coordinates we can assume that this subspace is V ect(e 1 , e 2 , . . . , e n ) thus there exist a vertex partition of Q n with set of translation say D. Then D ∪ {e n+1 , e n+2 . . . e n } define a vertex partition of Q n by translates of α(T ). 2

3 An example: antipodal paths

The antipodal vertex of a vertex x in Q n is the unique vertex x at distance n of

x. Notice that x=x + 1. An antipodal path is a path in Q n of n edges connecting some pair of antipodal vertices. We will say that an isometric embedding α(P ) of an antipodal path in Q n is canonical if v 0 = 0 and along the path the directions used are 1, 2, . . . , n in this order. We have thus, for any i ∈ {0, 1, . . . , n}, v i = i j=0 e j and θ W ( n i=1 λ i e i ) = ( n i=1 λ i i j=1 e j ). By a translation and a permutation of the coordinates we can always assume that an isometric embedding of P is canonical.

A vector u ∈ F n is of type 01, respectively of type 10, if there exists i 0 ∈ {0, 1, . . . , n} such that u = n i=i 0 +1 e i , respectively u = i 0 i=1 e i . Notice that 1 and 0 are the only vectors of both types.

2

We will now deduce theorem 8 from theorem 7 showing that this construction is in fact the analogue of Vasiliev's construction.

Alternative proof :

Let W = (v 1 , v 2 , . . . , v n ) and W = (v 1 , v 2 , . . . , v 2n+1 ),
where v i = i j=1 e j . By theorem 3 D n is obtained from a perfect code of length n Notice that, for any x ∈ F n and any a ∈ Z 2 we have:

C n = θ -1 W (d); d ∈ D n . Let λ : C n → Z 2 defined
θ W (x, 0, 0) = (θ W (x), 0, 0) + π(x).(0, 1, 1)

θ W (0, 0, x) = (0, 0, θ W (x))

θ W (0, a, 0) = a.(0, 1, 1).

Therefore:

θ W (θ -1 W (y), 0, 0) = (y, 0, 0) + π(θ -1 W (y)).(0, 1, 1),

θ W (θ -1 W (0, π(θ -1 W (y)) + λ(θ -1 W (d)), 0) = π(θ -1 W (y)) + γ(d) .(0, 1 , 1)
and θ W (0, 0, θ -1 W (y) + θ -1 W (d)) = (0, 0, y + d).

By linearity of θ W we obtain the expression of D 2n+1 .

2

Here also, if there exists u, v ∈ D n such that γ(u + v) = γ(u) + γ(v) then D 2n+1

is not a linear subspace. Notice that, if we set γ(u) = 0 for any u ∈ D n , we obtain a recursive construction of the linear subspace proposed by M.Ramras [START_REF] Ramras | Symmetric Vertex Partitions of Hypercubes by Isometric Trees[END_REF].

  defines a vertex partition of Q n by the embedded tree T . All vertex partitions of Q n by α(T ) can be obtained by this way.

Lemma 5 -

 5 If u and v are two distinct vectors both of type 10, or both of type 01, then u and v differ by a set of consecutive coordinates. If y + d = y + d then u + u = (y + y , 0, 0). But y + y = d + d and again by lemma 6 d = d , y = y and u = u . • γ(d) = γ(d ) thus u + u = (y + y , 1, d + y + d + y + 1). Therefore y + y is of type 01 and d + y + d + y + 1 must be of type 10, thus d + y + d + y of type 01. But d + d = (y + y ) + (d + y + d + y ) then by lemma 5 d and d must differ by a set of consecutive coordinates; a contradiction with lemma 6.

  by λ(c) = γ(θ W (c)).Consider C 2n+1 = (θ -1 W (y), π(θ -1 W (y)) + λ(θ -1 W (d)), θ -1 W (y) + θ -1 W (d)); y ∈ F n , d ∈ D n .By theorem 7 C 2n+1 is a perfect code of length 2n + 1 and thus by theorem 3 D 2n+1 defines a vertex partition of F 2n+1 by isometrically embedded antipodal paths whereD 2n+1 = {θ W (c); c ∈ C 2n+1 }.

Acknowledgment

The author is grateful to M. Kovše for pointing him towards the work of M.Ramras.

Proof : Assume first u = i 0 i=1 e i and v = j 0 i=1 e i for some i 0 , j 0 ∈ {0, 1, . . . , n} Assume w.l.o.g. that i 0 < j 0 we have thus v = u + j 0 i=i 0 +1 e i . For the second case notice that if u is a vector of type 01 then u + 1 is of type 10. Proof : Consider two translates of the embedded path say, x + α(P ) and y + α(P ).

If z is a common vertex of the two paths then for some i, j ∈ {0, 1, . . . , n} we have z = x + i k=1 e k and z = y + j k=1 e k . Thus by lemma 5 condition (ii) implies that the paths x + α(P ) and y + α(P ) are disjoint. By condition (i) every vertex belongs to some path.

2

If we assume that C is a linear subspace, the last condition is equivalent to the condition used by Ramras:

(ii') No element of C has as support a non-empty set of consecutive integers.

Let us recall Vasiliev's construction of perfect codes.

Theorem 7 (Vasiliev [START_REF] Vasilev | On ungrouped, close-packed codes (in Russian)[END_REF]) Let C n be a perfect code of length n. Assume 0 ∈ C n and let λ be a function from C n to Z 2 such that λ(0) = 0 and π be the parity function.

Then the set

Notice that if there exists u, v ∈ C n such that λ(u