Serge Monteix
email: serge.monteix@philips.com.

Yannick Le Maoult
email: lemaoult@enstimac.fr

Fabrice Schmidt
email: schmidt@enstimac.fr-arcens@enstimac.fr

Jean Paul Arcens

Quantitative infrared thermography applied to blow moulding process: measurement of a heat transfer coefficient

Keywords: blow moulding, IR camera, semi-transparent material, convection coefficient, anemometry

This paper deals with the heat conditioning stage of blow moulding process applied to P.E.T bottles forming. One of the most important stage of this process is the radiative heating stage which is performed with infrared ovens using powerful halogen lamps. To validate a 3D control volume thermal software, called Plastirad, developed in our laboratory, temperatures maps were needed on the plastic preforms as well as convective heat transfer coefficient inside the oven. This measurement has been performed with two different methods : IR thermography and hot wire anemometry. These two methods are investigated and the main results are compared to focus on the interest of IR thermography.

Introduction

The industrial injection stretch-blow moulding process of Polyethylene Terephtalate or P.E.T bottles (Rosato 89) requires a heating step prior forming. During this step, the amorphous preform corresponding to the form prefiguring the final bottle geometry is heated from the ambient temperature to the forming temperature between 105 to 115 °C. Due to the weak polymers thermal conductivity (0.25 W/m.K) and very high production rate (up to 40,000 bottles per hour); for the P.E.T, the heating step is made using infrared ovens [START_REF] Monteix | Modélisation du chauffage convecto radiatif de préforme thermoplastique pour la réalisation de corps creux[END_REF]. These oven are constituted with a row of special halogen-lamps, which permit a rapid heating with high flux intensities (300-400 kW/m 2 ) and a cooling fan to control the external surface temperature and to avoid the thermal crystallisation which would make the transparent amorphous preform completely opaque and more fragile. Then, the ventilation appears as a key parameter controlling the preform temperature distribution. Our investigation deals with the measurement of this heat transfer coefficient using the experimental set-up we reproduced from the industrial machine to the laboratory scale. To do this, two methods were used: one based on measurements made on the preform surface with an IR camera, the other using conventional anemometer technology to measure the air velocity and temperature. Then comparison of these two methods is analysed and discussed.

Experimental set-up

Figure 1 shows the experimental set-up developed in our laboratory, constituted with six halogen lamps (Philips 1 kW-235V) and a flat aluminium reflector located behind the lamps. The cooling fan disposed opposite to the lamps is enclosed inside an aluminium box and the flow produced by the fan passes through a set of slits and is blown onto the preform surface placed normally to the air jet. Due to the smooth aspect of the aluminium box, we consider this part of the set-up as a secondary infrared reflector. The preform, 100 mm height and 25 mm in diameter (Fig. 1), can rotate axially in the middle of the oven, this point is of prime importance for a uniform heating through the preform circumference. However the narrowness of this experimental set-up leads to instrumentation difficulties for the different sensing probes.

Heat transfer coefficient computation using an infrared camera.

Parameters influencing the measurement accuracy using an infrared camera and the P.E.T preform.

The heat flux received by the camera is related to the radiative properties of the materials observed; here, the P.E.T preform but any elements in the thermal scene (Fig. 2) such as reflectors, lamps can contribute to the total heat flux integrated by the camera. Thus, precautions are needed in order to measure the real temperature of the preform with a good accuracy. So we had to describe the environment participation; this can be achieved, at first, by the radiometric equation .

) ( . ) ( . ). ( ) ( 0 0 0 env pref pref pref atm atm eq T V T V T T V [1]
In equation [1], V 0 (T eq ) corresponds to the camera output voltage; T eq , the temperature corresponding to an equivalent blackbody temperature obtained with the system calibration; atm , the atmosphere transmission coefficient between the preform and the camera. This coefficient is assumed to be equal to one in our application due to the very short distance between the camera and the preform : about 0.5 m; pref and pref are respectively the emissivity and reflectivity of the P.E.T preform integrated over the spectral band of the camera chosen for the measurement. This implies a study of these two parameters to choose the more relevant spectral band. T pref is the real preform temperature and T env , the equivalent blackbody environment temperature. P.E.T spectral properties measurements were obtained using 3-mm thick sheets. These sheets were obtained by injection moulding in the same conditions than tubular preforms. Differential scanning calorimetry tests (DSC) were carried out in order to measure the sheets crystallinity. Samples of 10-15 mg weight were cut away from the thick sheet. The DSC tests were performed on a Mettler calorimeter at a heating rate of 10°C/min. The average crystallinity value of these measurements was about 5 %, corresponding to an amorphous P.E.T as well as experiments performed on tubular preforms, which confirms that diffusion of the radiation can be neglected. Then the P.E.T spectral properties measurements have been made using a Perkin-Elmer FT-IR spectrometer for the range of wavelengths in [1 -25] µm. The initial 3mm thick sheets were polished in order to obtain samples of 0.3 mm thickness. We assumed in agreement with (Weinand 87) that these properties are independent of the temperature variation. From apparent transmission measurements Ta for the range of thickness [0.3-3 mm] and (Howell 92), we get : , i.e These measurements lead to an intrinsic absorption coefficient K a independent of d (Beer's law is relevant) and reveal that P.E.T presents successively transparent and opaque region over large spectra between 1 to 25 m (Lebaudy,89), [START_REF] Monteix | Comparison between a numerical model and an experimental approach of preform infrared heating -recent results[END_REF]. We have then computed the Planck's mean absorption coefficient with [3] :

i d T I i d T I K K T K K b K b a a ) , ( ) , ( .
) 400 ( 400 , 400 ,

[3] T = 400K is the classical processing temperature of the perform and I b the spectral distribution of blackbody emission. Equation [3] gives K a (T) = 200 m -1 for i = [1-3 µm] region, 4000 m -1 in the [3-5µm] region and 30000 m -1 in the [5-20µm] region. Using the optical thickness K a .d, and K a .d = 1, we got a photon mean free path d = 1/K a to be compared to the thickness of the preform (3 mm) : respectively d = 5 mm for the [1-3 µm] region, d =0.25 mm for the [3-5µm] region and 30 µm for the [5-20 µm] region. This last value corresponds to an opaque material for which radiation come only from the material surface. It means correction used by [START_REF] Hajji | Radiation pyrometry on semitransparent sheets I:Gray Media[END_REF] for semitransparent polymers are not required here and measurements using a Long Wave spectral band [8-12 µm] are surface temperature measurements. For this reason, we have chosen an AGEMA 880 Long Wave infrared camera connected to a real time 12 bits thermography software among the three systems we get in Ecole des Mines: Very Short Wave (VSW) [0.75-1.1 µm], Short Wave [3-5 µm] (SW) and Long Wave (LW). The corresponding values of emissivity and reflectivity calculated for the radiometric equation [1] in order to compute the real preform temperature are respectively = 0.93 and = 0.07. Then the camera was calibrated versus the temperature using a Land blackbody between 30°C and 175°C corresponding to the full window of the phenomenon we wanted to observe. Ten points over this range were recorded and led to an accuracy better than 0.5°C in comparison with the blackbody reference temperature.

Oven environment influence on the IR measurements: a first approach.

To control the temperature measurement of the preform inside the oven, we have firstly estimated the different irradiances incoming to the lens of the camera from the different areas of the oven, mainly : E1 : irradiance received by the camera from the preform E2 : irradiance received by the camera from the quartz tube and reflected by the preform E3 : irradiance received from the ventilation system and reflected by the preform.

Assuming that all surfaces are diffuse for reflection and emission, composed of grey materials, the different terms, neglecting the atmosphere absorption and multiple reflections are given by :

c p c p v v v v v p c v c p c p Q Q Q Q Q p c s c p c p p p p c pref F S F T f T S S E F S F T f T S S plamp E F S T f T S S E ) ) , ( (
) ) , ( ( ) , ( 4 3 4 2 4 1 [4]
and E t =E 1 +E 2 +E 3 is the total irradiance seen by the camera. f( ,Tp) corresponds to the fraction of radiation emitted between the two wavelengths 1 and 2 , here: 8 and 12 µm. These expressions depend also on geometric interaction depicted by view factors computation F; S c is the input surface of camera lens, S Q is an equivalent area corresponding to the halogen lamps, S V is the effective area of the ventilation support, S P the effective area of the preform and F i-j are the different view factors between lamps, preform and ventilation support (example: F Q-P is the view factor between lamps and the preform). In each case, to compute these view factor using analytical formulation [START_REF] Siegel | Thermal Radiation Heat Transfer[END_REF], the real geometry is replaced by an elementary one (figure 2c). These assumptions lead to the computation of simple view factors between two rectangles with dimensions reported in table 1.

Then, with [4], we compared the relative contributions of each term to the total irradiance. The irradiance of each source were computed versus its own temperature (table 2) considering maximum value at the cooling system restart. These comparisons lead to conclude that the oven environment contributions are negligible. Using [4], we have simulated, on figure 3, the equivalent blackbody temperatures in the different cases : superimposing the blackbody reference and the irradiance received by the camera with all the perturbation (E t = E 1 +E 2 +E 3 ) without any correction (assuming the preform is a black body) and then taking into account the optical P.E.T properties and an environment temperature equal to 25°C. This computation leads to a discrepancy of 1°C for a 120°C preform temperature reference (figure 3).

Geometry

This numerical approach has been compared to experimental data; a sequence of images has been recorded when the preform was just entering inside the oven : it was possible to measure a new equivalent blackbody temperature related to the cumulated reflections on the preform just before heating. This temperature was close to 28°C. The difference between the new and default environment temperature (25°C) is 3°C, slightly different from the computed value (1°C), but we assumed that a correction was not necessary.

Measurement of the heat transfer coefficient-methodology

The computation of the heat transfer coefficient distribution along the preform height is based on the IR images acquisition (Figure 4).

The acquisition rate was equal to one image per second for a preform in rotation inside the oven. More particularly we dealt with the spatial average temporal temperature variation over 7 circular zones regularly disposed along the preform height. Here, each circular zone had a diameter of 5 mm, this size is a compromise between an area large enough to check the homogeneity of temperatures and spatial resolution of our camera. Then three steps are required for the measurement : -The first step : a 20s heating stage which permits to get a sufficient preform temperature; the ventilation system is on to avoid the crystallisation of the perform surface.

Figure 4. Areas of interest and methodology

-The second step : an equilibrating phase of 10 s is applied to leave the preform temperature evolve freely due to conduction , radiation and natural convection in order to tend to a zero gradient through its thickness.

-The third step : the ventilation is switched on. We analysed more particularly this phase.

The 10s equilibrating step we have imposed to the preform was shown to be appropriate by measurements of the inner and outer preform temperature [START_REF] Monteix | Comparison between a numerical model and an experimental approach of preform infrared heating -recent results[END_REF].

Applying these conditions, we assumed that, at the ventilation restart, the preform presented for each measurement target along its height, a uniform temperature through the thickness corresponding to the surface temperature at this time. In this scope, we have chosen to analyse the temperature distribution using a transient one-dimensional semi-infinite analytical model (DeVriendt, 1984) through the preform thickness by solving the following equation :

t T x T 1 2 2 [5]
with the following initial and boundary conditions :

ini 0 T t) , (x ; )) , 0 ( ( T ; ) 0 , ( T t T T h x T t x T g x ini [6]
T ini is the initial temperature of the perform at the beginning of the cooling and T the temperature of the air far away from the perform. Finally, solving this equation leads to the temperature relation ship for temperature versus time t, thickness x and h :

)] 2 ( . [ ) 2 ( ) , ( 2 2 t h t x erfc e t x erfc T T T t x T g t h x h ini ini g g [7]
where is the thermal diffusivity of the material and here, its thermal conductivity. For P.E.T these values are = 1.4.10 -7 m 2 .s -1 and = 0.25 W/m.K. In this case, the heat transfer coefficient h g is the amount of the radiative and convective transfer. The convective part is obtained after subtracting the radiative one and is given by :

env env r g c T T T T h h h ) ( - h 4 4 g , [8]
where T env is a radiative environment temperature.

One-dimensional analytical model accuracy

Since the preform geometry corresponds to a cylinder, the approach used is a rough assumption and to evaluate the resulting discrepancy, we have compared temperature computations from this one-dimensional model to ones issued from a 3D Control Volume Model (CVM) using "Plastirad" software dealing in this last case with a tubular geometry of 3 mm thickness close to the preform one. More particularly, we compared temperature evolution, during a cooling step, from a uniform initial temperature, due to three different heat transfer coefficients of 5, 35 and 65 W m -2 K -1 . We have also taken into account the same P.E.T properties and a constant environment temperature T env equal to 30 o C. We have superimposed on figure 5 the external surface temperature T ext and the internal one T int at a 3 mm depth versus the time for the greatest value of h g .

We can notice that, for the internal temperature, a difference is noticeable between CVM and 1D computation for a time close to ~7s. This time corresponds, for the thickness e equal to 3 mm, when the surface perturbation is viewed at 1% on the internal surface of the preform, which leads to (DeVriendt 1984), t ~(e 2 / 13.24 ~7s. This divergence is well described by the CVM. To quantify the discrepancy between the two approaches, we have computed the relative difference for each time step and from these value we have computed the average relative difference versus the time for each h (Table 5).

Figure 5 . 1D temperatures computation and CVM

; h g = 65 W m -2 K -1 hg (W m -2 K -1 ) Error / T_int (%)
Error / T_ext (%) 5

1.1 0.7 35 2.9 1.8 65 2.8 4.4

Table 5 : Analytical model accuracy Differences between temperatures computed respectively with the CVM and the 1D analytical method are low: the mean relative error is less than 2.3%. These value lead to conclude to the interest of the simple 1D semi infinite approach which gives a good description of the preform surface temperature variation versus the heat transfer coefficient for times close to 7 s.

Heat transfer coefficient computation

Now, using the 1-D approach and surface temperature values extracted from each interest zones, we processed to the heat transfer coefficient identification minimising the total relative difference between experimental and temperature computations, assuming a constant environment temperature. We have superimposed (figure 6) the measured and computed variation of the surface temperature for each zone along the preform height. These curves correspond to an air flow temperature equal to 20 o C at the beginning of the cooling step, after 20 s heating and 10 s of equilibrating phase. This air flow temperature was measured with shielded iron-constantan thermocouples as shown on figure 1 (two thermocouples are visible near the perform); compared to the effective duration (for the computation of h c ) of the cooling phase, no significant variation was observed on this temperature, in fact equal to the ambient temperature of the laboratory during the experiments. The computation of the heat transfer coefficient was, then, done for a cooling time of 7 s according to the trends viewed on the figure 5 (for this specific time the 1D model is valid). The difference between measured and computed values is very low : minor than 2 °C, around 110 °C.

Figure 6. Comparison between measured and computed preform temperatures

The sensibility to the air flow temperature is plotted below : variations of the heat transfer value corresponding for each air temperature. On figure 7, the convective heat transfer coefficient is plotted. We saw that variations of the air flow temperature tend to present high variations in the computation of the heat transfer coefficient. On figure 7, 20°C is a measured temperature and due to the difficulty to get accurate measurements on this air temperature (intrusive effect of the temperature probe in the flow, recirculation zones), only the effect of a presumed variation around 20 °C was reported for observed air temperatures in industrial workshops: 15 °C to ~ 40°C. Then to keep the same cooling profile than for an air temperature of 20 °C (Fig. 6), the new heat transfer coefficient for an air temperature of 40 °C must increase (by analogy with the Newton law of cooling: q = hc(Tsurface-Tair), for q assumed to be a constant and Tair increased, h is increasing too); in this case, that can be possible only by increasing the flow velocity as it can be shown with classical convection correlations. It reveals here the importance of the air flow temperature input data in order to increase the accuracy of the heat transfer computation.

Figure 7. h c computation versus the air temperature value (sensibility analysis)

Heat transfer computation from air flow velocity measurement

This second method is based on the air flow velocity measurement near the preform. Due to the narrowness of the space between the lamps and the oven (less than 10 cm), an ANS Snelco hot wire anemometer with a speed range of 0-30 m/s was used at different points along the preform height and for a same level two locations were used as shown on the figure 8, the distance between the sensing part of the anemometer and the preform was equal to 10 mm for the two locations and for the front position, 35 mm from the outlet slits respectively. For each location, several measurements of the air speed were made; at last the average value of the speed for all the measurements of a same level was computed due to the fact that the Plastirad software can only take into account an average hc coefficient versus the height of the preform.

Measurements were performed lamps switched off at the ambient temperature (20°C) because of the high radiative heat flux from the halogen lamps during the heating stage.

Using the air velocity measurement, we computed the heat transfer convective coefficient h c with the Nusselt number, taking into account air thermophysical properties at the ambient temperature (~20°C).

Computation of the convective heat transfer coefficient

The Reynolds number Re = V.d p / , where V is the air velocity , d p the characteristic length of the air flow : here, the preform diameter (d p =25.5 mm) and the air Kinematic viscosity. In our case, the air velocity varied between 1.5 to 5 m.s -1 (the lowest part of the scale of our anemometer). Then the Reynolds number was defined by the following interval: 4000 < Re < 10000, so we assumed that the air flow was turbulent and used, in a first approach (because in this particular situation, the flow was not a uniform flow), the Hilpert correlation from the literature [START_REF] Incropera | Fundamentals of Heat and Mass Transfer[END_REF] giving the Nusselt number versus Reynolds and Prandtl corresponding to this particular configuration : [9]

Finally the air velocity measurements and the corresponding heat transfer distribution along the preform height are plotted on figure 9.

These results show a characteristic variation due to the specific ventilation system we used and a high measurement dispersion (in accordance with metrology advises, the uncertainty bars for the computation of h c are twice the standard deviation of the same parameter): the uncertainty for the h c coefficient is equal to 29 % for the maximum value and 23% for the minimum one; but we want to mention here that it was particularly difficult, inside the oven environment, to visualize air flow using smokes or other techniques due to the halogen lamps and electrical wire proximity. To compare the two methods, we have superimposed results of the h c computations. Results obtained by the thermography method represent the h c average value obtained for different presumed air temperatures resulting of the sensibility analysis around air temperature of 20°C . The corresponding uncertainty bars results also from twice the standard deviation of these values (figure 10).

We note, here, that the two methods lead to a same heat transfer coefficient scale value. However the distribution along the preform height is different. We note a large uncertainty of the two methods due to the high non uniformity of the air flow and the temperature variation inside the air flow. Also, the 1D method tends to integrate and to smooth the values of h c .

Figure 10. Comparison of h c along the preform height obtained by the two methods

Conclusion

In this paper, we developed different methods to estimate the heat transfer coefficient along the preform height coupled with an air cooling system. At first, we analysed the environment participation inside the oven and showed that, despite all the reflectors inside the oven, the participation of the preform environment was negligible in a first approach. Then, we developed a specific analysis involving IR thermography using images recorded all along the preform height and focused this analysis on the cooling step in order to simplify the problem. A 1D analytical model was used to compute the preform surface temperature distribution versus the heat transfer coefficient. In order to validate this computation, we compared the accuracy of the 1D model with a CVM. This comparison leads to conclude that temperature discrepancies are less than 10%. Finally, for each zone of the preform, we measured a local and average heat transfer coefficient h c (z). The comparison between the thermographic measurements and values computed with a classical anemometric method provided heat transfer coefficients in the same range. However, we noted that the two methods presented a large uncertainty. The method based on the thermography seems very sensitive to the air temperature and takes into account all the components of the air flow that could explain these discrepancies. So, due to the relative proximity of the heat transfer coefficient obtained by these methods, we will emphasize on the interest of the thermographic method which deals with all the preform height and limits the measurement time. These measurements were also tested successfully on industrial machines and led to realistic values of h c which were in the range of 35-60 W/m 2 .K, range very close to our laboratory results but smaller than the values found by [START_REF] Lebaudy | Etude et simulation de la répartition des températures dans un cylindre creux de P.E.T soumis à un rayonnement infrarouge[END_REF]: ~ 80 W/m 2 .K. These values were obtained with an iterative procedure and a 2D axisymmetric finite difference method. In order to confirm potentiality of our method, it would be interesting to reproduce this comparison dealing with a bigger heat transfer coefficient range to test the effect of the variation of h c on the quality of this specific process. Also, to increase the accuracy of the thermographic method, a particular attention will be attached to the air flow temperature measurement enhanced by the radiative environment and intrusive action of the thermocouples. In the same way, if a real reference method, based on a more realistic correlation for jets, is needed to "calibrate" the thermographic method, a work has to de done to provide a better description of the flow inside the oven; for that, a CFD software such as FLUENT could improve our knowledge of the ventilation system. At last, we are currently working on a specific ray tracing method devoted to accurate radiative transfer computation inside the infrared and also to optimise the position of the IR camera to be less sensitive to reflections.

Figure 1 .

 1 Figure 1. Experimental set-up and infrared oven

Figure 2 .

 2 Figure 2. The thermal scene viewed from the infrared camera

Figure 3 .

 3 Figure 3. Equivalent blackbody temperatures.

Figure 8 .

 8 Figure 8. Anemometer set-up description. Left: upper view, right: front view

Figure 9 .

 9 Figure 9. Air speed velocity measurement and computation of h c

Table 1 .

 1 Dimension of the thermal scene for analytical computation of view factors

		W		H	Distance/preform	View factor
		(mm)	(mm)	(mm)
	Preform	5		100	0
	6 halogen lamps system 122.5	100	34	FQ-p=0.989.10 -3
	Cooling system	86.5		100	34	Fv-p=0.127. 10 -2
	Radiative source	T (K)			Comment
	Preform	373	Average temperature measurement; after the cor-
		responding heating step and during natural cooling
	6 halogen lamps :	730 Current temperature computed during halogen lamps
	quartz temperature		running. Certainly higher than during the cooling
					stage
	Cooling system	320	From thermocouples stickled over the ventilation
				support and shielded inside aluminium cavity

Table 2 .

 2 Temperature values of direct and indirect radiative sources

	Radiative	Properties [8-12] m	Source of Data
	source
	Preform	= 0.93, =0.07	measurements
	the 6 lamps	=0.82	(Monteix, 2001)
	Cooling sytem	=0.52	Measurement at 45 o C using the camera and a black body at 100 o C

Table 3 .

 3 Radiative properties of thermal sources

	Radiative source	i W/m 2	µm	Ei/Et (%)
	Preform	0.019		92.2% /direct
	6 halogen lamps system	0.0016		7.6%, reflection
	Cooling system	4.03 E-5		0.2%, reflection

Table 4 .

 4 Estimation

of the different irradiances received by the camera

The optical properties in

[8][9][10][11][12]

µm are reported in Table

3

and comparisons in table 4.

Acknowledgements

The authors would like to thank The Nestlé Waters Company for participation in this work, and specially G. Denis for his technical support.


F Deschaux-Beaume

F Schmidt

N Fre Âty

T Cutard

J C Boyer

C Levaillant

Numerical simulation of failure prediction for ceramic tools: comparison with forging test results

Keywords: Finite-element method, Fracture prediction, Ceramics, Forging

A fracture prediction criterion for brittle materials has been introduced in the POLLUX ®nite-element code in order to predict the risk-ofrupture of ceramic tools during a forging operation. POLLUX is a software dedicated to the simulation of forging operations, initially developed by INSA (Lyon). The chosen probabilistic fracture model is based on the weakest-link theory and the statistical theory of Weibull. A surface approach or a volume approach can be retained on the basis of the type of critical ¯aws in the ceramic. Two different criteria are available in order to characterise the stress state, considering the tensile normal stresses and neglecting the compressive stresses. An identi®cation procedure of the critical ¯aw type is presented for a particular ceramic material. Statistical parameters of ceramic fracture have been determined experimentally using bending tests performed under environmental conditions close to those of forging. A constitutive equation of the workpiece material has been proposed, issued from torsion tests. In order to validate the model in the case of ceramic tools subjected to multi-axial stress states, a particular con®guration has been de®ned to compare the simulation predictions with the experimental results. A forging test has then been developed, in which a billet of superalloy is formed in a ceramic tool up to its fracture at the temperature of 1423 K. The experimental distribution of tool fracture, according to the strain of the billet, is in good agreement with fracture predictions computed by the simulation.

Introduction

Isothermal forging is a very effective forming process for some refractory alloys, such as nickel-based materials. This technique should allow the production of near-net-shape parts directly from preforms of simple shape. However, the process is based on the superplastic properties of the forming materials, and it must be carried out at high temperature and low strain rate. This means that the material for the forging dies, which will be used at the same temperature as the workpiece, must conserve high mechanical properties at high temperature. The forging tools are actually made in Mo-based alloys, but the high reactivity of this element with oxygen imposes the carrying out of the process under an inert gas. The use of ceramic dies should be a less expensive way as the good oxidation resistance of these materials does not require a controlled atmosphere.

However, ceramics are inherently brittle. This characteristic leads, in particular, to a wide variation in the material strength. The usual design techniques are, therefore, not suited to ceramic components, and statistics and reliability analysis are needed. Statistical fracture models based on the weakest-link theory have been developed previously. The ®rst probabilistic approach used to account for the scattering in the fracture strength of brittle materials was introduced by Weibull [1]. This model was proposed initially for constant uniaxial stress ®elds. New criteria were added later to extend the Weibull approach to multi-axial stress ®elds [START_REF]Fracture[END_REF][START_REF] Evans | [END_REF]. Another model, developed by Batdorf, combines the weakest-link theory and linear elastic fracture [4].

Except in some simple cases, it is necessary to use numerical simulation in order to predict the risk of failure of ceramic parts with these models. Some studies have been conducted previously in which brittle fracture models were coupled to ®nite-element simulations in order to calculate the failure probability of particular ceramic structures. In some cases, a post processor has been developed for a speci®c software to analyse the probabilistic behaviour of simple ceramic structures [START_REF] Berdin | Proc. of Fracture Processes in Brittle Materials[END_REF] or of uniaxially-loaded ceramic composites [START_REF] Guillaumat | Proc. of Journe Âe nationales sur les composites[END_REF]. However, most of the studies in this ®eld have been made by the National Aeronautics and Space Administration [7±9]. A computerised design program for the reliability evaluation of ceramic components has been developed from this research work. This software can be coupled to the ANSYS, ABAQUS, NASTRAN and COSMOS commercial ®nite-element programs. Nevertheless, all of these programs are inconvenient for the simulation of the forging process as large strains take place in the workpiece.

The main objective of this study is to develop a ®niteelement model for the failure prediction of ceramic forging dies as a function of the service conditions in order to improve the design of these components.

Description of the failure prediction software

Brittle fracture modelling

The wide variation in the strength measured for brittle materials is due to the nature and distribution of intrinsic microscopic ¯aws. These ¯aws are unavoidably present as a result of the material manufacturing and machining operations, or of a damage process. Distinction must be made between volume and surface ¯aws, which lead to different strength distributions. In the case of a coarse material microstructure, the volume ¯aws are generally more critical than the surface ¯aws. However, the surface ¯aws may be more critical, for example, in the case of rough machining or of oxidised surface states. The identi®cation of the nature of critical ¯aws allows one to choose the approach to the analysis. If a surface ¯aw is chosen, the volume ¯aws are assumed to be not responsible for the fracture.

As the ¯aw distribution cannot be determined, some statistical functions are used to describe the strength variability of a component. The Weibull model has been chosen in the present study because the statistical parameters are easy to determine experimentally. The Weibull theory associates a strength c to each elementary volume Ve (or elementary surface Se for the surface approach). The strength of the material is characterised by the distribution function of c , denoted F(), which is de®ned as the number of all elementary volumes (surfaces) having c , divided by the total volume (external surface). This function also gives the probability P of randomly choosing an elementary volume (surface) having a strength equal to or less than , or also the probability that this volume fails under the stress :

P c F (1)
Any distribution function may be written in the form

F 1 À expÀ (2)
where () must be a positive and non-decreasing function, vanishing at a value u , called the threshold stress. The advantage of this equation appears when considering the failure probability of the structure constituted by all the elementary volumes (surfaces), using the weakest-link theory. When supposing that the total volume V (surface S) fails if any one of its parts fails, then the probability of non-failure P NF , under a uniform uniaxial load , is equal to

P NF 1 À P n expÀn ( 3
)
where n is the number of elementary volumes (surfaces) of the structure (V nVe or S nSe).

The failure probability P n of the structure takes then the simple form

P n 1 À expÀn (4)
Weibull chose a power form for the function (), which is the simplest function satisfying the general conditions described previously:

À u 0 m (5)
where u , 0 and m are parameters depending on the characteristics of the material.

The threshold stress u is usually taken as equal to zero. The two other parameters can then be identi®ed easily by comparison with the results of fracture tests.

In the case of a non-uniform uniaxial load , the failure probability takes the integral form

P v 1 À exp À 1 Ve Z V 0 m dV 0 @ 1 A (6)
for the volume approach, or

P s 1 À exp À 1 Se Z S 0 m dS 0 @ 1 A (7)
for the surface approach.

The previous relationship can be generalised to multiaxial stress states. Two strategies are available: the ®rst consists of considering an equivalent stress that characterises the loading at each point, this criteria sometimes leading to unsafe estimates of the failure probability; whilst the second consists of assuming the independence of the actions of each principal tensile stress [START_REF]Fracture[END_REF].

Taking into account the ®rst approach, Weibull proposed the calculation of an equivalent stress by averaging the normal tensile stresses in all directions. This model leads to a complex integral formulation, which may increase the computation time.

In order to obtain shorter computation times, only the maximal normal stress has been considered in the present model, as has been done by other investigators [START_REF] Berdin | Proc. of Fracture Processes in Brittle Materials[END_REF]. This variable is also the maximal principal stress at each point of the structure.

In the second approach, the fracture probability depends on all the tensile principal stresses. The independence of the actions of the principal stresses means that the reliability of the structure is the product of the reliability of the individual principal stresses acting independently. The fracture probability is then given by the following relationship:

P v 1 À exp À 1 Ve Z V 1 0 m 2 0 m 3 0 m dV 0 @ 1 A (8)
for the volume approach, or

P s 1 À exp À 1 Se Z S 1 0 m 2 0 m 3 0 m dS 0 @ 1 A (9)
for the surface approach, where 1 , 2 and 3 are the positive principal stresses or take the value zero for the compressive principal stresses.

It is interesting to note that in both models, the compressive stresses are not taken into account. The fracture is supposed to be always the result of tensile stresses. This choice has been made because the compressive fracture strength of ceramics is three or four times greater than the tensile strength. It is very unlikely that a structure fails under the effect of compressive stresses.

Numerical introduction of the model

The model described previously has been introduced as a subroutine in the POLLUX ®nite-element code [10]. POLLUX is a specialised software developed by INSA-Lyon [11] in order to predict the stress state in the workpiece and the surrounding tools during a forging process. First, the thermo-mechanical analysis is achieved, leading to the calculation of the stress tensor at each point of the tools, which are supposed to show a perfectly elastic behaviour. Then, the failure probabilities of the tools are determined, using the Weibull parameters and the criteria chosen by the user.

If Ve or Se are arbitrarily taken to be equal to 1, the discrete decomposition of the models gives

P v 1 À exp À X iel X ig 1 0 m Vig ! (10)
for the volume approach, or

P s 1 À exp À X inc 1 0 m Sinc ! (11)
for the surface approach if the maximal principal stress criteria (MPS) is chosen, and

P v 1 À exp À X inc X ig X 3 i1 1 0 m Vig ! (12)
for the volume approach, or

P s 1 À exp À X inc X 3 i1 i 0 m Sinc ! (13)
for the surface approach if the principle of the independent actions (PIA) is chosen.

1 , 2 , 3 are the principal stresses in the Gauss point ig of element iel or in the contour node inc ( 1 ! 2 ! 3 ) or are equal to zero for the compressive principal stress. V(ig) is the volume associated with the Gauss point ig, and S(inc) is the surface associated with the contour node inc.

The elementary failure probability DPr, or local risk-ofrupture, is also calculated at each Gauss point in the case of the volume approach, and interpolated to the nodes. This variable corresponds to the fracture probability of an elementary volume of constant size equal to 1, subjected to the stress state of the point considered. The values can be stored in a data ®le for a graphical display of the critical regions. Another data ®le, containing the failure probability of the overall structure as a function of the time or tool displacement, is also generated. At the end of each increment, the program compares the failure probability of the tools to the maximal permissible probability chosen by the user. The computation is interrupted when this limit or the maximal tool course is achieved.

Reliability analysis procedure for ceramic forging tools

Identification of critical flaw type

Reliability analysis of ceramic structures with the software described previously ®rst requires a failure criteria to be chosen corresponding to the material and the use conditions. If the choice between MPS and PIA criteria cannot be made by considering physical observations, the choice of the type of approach (surface or volume) must be based on these observations.

In the case of a forging process with ceramic tools, the material strength is modi®ed because of the damage inherent to the use conditions. Two types of ¯aws can be created: surface ¯aws, essentially produced by the oxidation of the material at the forging temperature; and internal ¯aws (microcracks), which can be formed because of the mechanical loading of the tool or because of the thermal cycles.

In order to identify the critical ¯aw type formed during the forging process, the fracture strength of the initial material and the material subjected to a damaging factor (oxidation, mechanical or thermal cycles) are compared. The failure probability of the material, as a function of the stress, is determined using four-point bending tests (see Section 3.2) at room temperature. Fig. 1 shows the experimental points and the corresponding interpolated Weibull functions for the different states of the material. First, it is to be noted that the oxidised material has a lower strength than the initial material. Then it appears that mechanical and thermal cycles have no signi®cant effects on the material characteristics. The strength distribution of the material after mechanical cycles is, indeed, very close to that of the initial material, and the strength distribution of the material after oxidation and thermal cycles is approximately the same as for the uncycled oxidised material.

This means that the surface ¯aws, created by the oxidation of the material, are more critical than the initial ¯aws, or the internal ¯aws resulting from cyclic loading. According to this result, the surface approach must be chosen for the calculation of the reliability of the ceramic tools during the forging process.

Determination of statistical parameters

The second step of the reliability analysis of the ceramic tool consists of the experimental determination of the statistical parameters of the Weibull model corresponding to the material and the processing conditions. Four-point bending tests (Fig. 2) are carried out at 1423 K on a set of N specimens oxidised previously at 1423 K for 50 h. The failure probability of each sample is given by the equation

P ex i N 1 ( 14
)
where i is the rank of the maximal stress max in the corresponding specimen at the fracture point. max is given by

max 3 2 FL À 1 bh 2 (15)
The Weibull distribution corresponding to the bending tests, in the case of surface critical ¯aws, can be expressed as follows:

P th 1 À exp ÀS eff max u m ( 16
)
where S eff is the effective surface determined from the surface in tension, S, of a specimen:

S eff Z S max m dS ( 17
)
Eq. ( 16) can also be written as

ln ln 1 1 À P th ln S eff m u m ln max ( 18
)
The Weibull plot, which expresses the evolution of ln(ln(1/(1 À P))) as a function of ln( max ), is usually used to represent the experimental results. The Weibull parameters are estimated by ®tting Eq. ( 18) in the experimental points of a Weibull plot, using the least-squares method (Fig. 3). The Weibull parameters m and u obtained are equal to 12.4 and 580 MPa, respectively.

Determination of the workpiece rheology

A further set of data required for the reliability analysis is the rheological behaviour of the workpiece material. Under forging conditions, the superalloys have a superplastic behaviour. In consequence, the ¯ow stress depends on the strain rate. In order to determine the constitutive equation of this material, torsion tests have been carried out at the forging temperature with three different strain rates ranging from 0.03 to 0.3 s À1 . During the test, the experimental torque is measured as a function of the turn number T. The equivalent strain rate " " and strain " ", calculated at the periphery of the specimen, are given by "

" 2 TR 3 p b with T dT dt (19) " " 2 TR 3 p b (20)
where R and b are the radius and the length, respectively, of the specimen. The Von Mises equivalent stress " is given by

" 3 p 2R 3 ñ m (21)
where ñ and m represent the sensitivity of the torque to the strain and to the strain rate. The evolution of " with " " is plotted in Fig. 4 for the different strain rates. The curves show a stress peak (" p ), followed by softening and stabilisation of the stress values (" s ).

According to Berdin et al. [12], the stress peak observed in this type of material at low strain levels is due to a dynamic recrystallisation phenomenon. The material behaviour can be approximated by the model

" 1 À " p " s (22)
with

" p A p " " mp ; " s A s " " ms ( 23
)
where the different constants are determined by ®tting the model to the experimental results. mp and ms are the strain rate sensitivity coef®cients, equal to 0.22 and 0.15, respectively, whilst A p and A s depend only on the temperature and are equal to 150.1 and 109.6, respectively. represents the recrystallized fraction, given by an Avrami type equation [12]. Fig. 4 shows a good agreement between the experimental results and the curves calculated from this model.

4.

Confrontation with the forging-fracture test results

Description of the forging test

In order to validate the model and the approach described previously, fracture tests have been carried out on ceramic tools under near-forging conditions, the results being compared to the numerical results of the tool reliability. Fig. 5 shows the samples chosen for this confrontation. The ceramic and the workpiece are heated to 1423 K in an induction furnace. During the test under isothermal conditions, a superalloy workpiece is pressed by a punch into the ceramic specimen until its fracture. The forging load and the punch displacement are measured throughout the tests.

After the tests, realised under the same conditions on a set of ceramic specimens, the failure probability calculated from Eq. ( 14) is expressed as a function of the punch displacement. A program of 10 fracture tests has been planned to obtain an experimental distribution of ceramic fracture.

Failure probability calculation

Good agreement between the experimental results and the numerical simulation predictions has been observed. First, the calculation of the stress state in the specimen has shown that the orthoradial stress is the maximal component of the stress tensor in the ceramic (Fig. 6). This predicts that, at fracture, the cracks in the ceramic specimen must be longitudinal: this has been veri®ed experimentally and can be observed clearly in Fig. 5. Secondly, the failure distribution of the ceramic specimens, expressed as a function of the punch displacement, obtained from the numerical calculation are presented in Fig. 7. The experimental values measured during the forging tests are very close to the calculated values. Both MPS and PIA criteria available for the calculation have been tested in order to determine as to which is more ef®cient for the present case. Nevertheless, the choice seems to be dif®cult because of the weak difference obtained from the two criteria. This is so because, for the ceramic geometry chosen, one of the principal stresses is always predominant, and the other two have no signi®cant effects on the risk-ofrupture. Thus, except for quasi-equi-biaxial stress states, both the MPS and the PIA criteria are ef®cient in predicting the failure probability of the ceramic material chosen.

Conclusions

A fracture-prediction program for brittle materials has been introduced in the ®nite-element code POLLUX in order to predict the risk-of-rupture of ceramic tools during high temperature forging operations. For ef®cient use, the following characteristics must be known: (i) the type of the critical ¯aws in the ceramic; (ii) the statistical parameters of the fracture strength distribution of the ceramic; and (iii) the metal workpiece rheology. These characteristics can be measured by classical mechanical tests (bending, torsion or tension tests).

The software is then a very useful prediction tool for the fracture of ceramic parts as: (i) it permits the prediction of the failure probability of a ceramic tool with a given geometry, as a function of its use conditions; and (ii) it can also give elements to predict the localisation and the orientation of the fracture cracks. So, such a software can be very helpful in the design of the geometry and of the dimension of large ceramic tools.

The comparison between the numerical simulation results and the forging fracture tests has shown a good agreement. Nevertheless, the test con®guration has not permitted differentiation between the two available criteria (MPS and PIA). Complementary experimental tests could be performed in order to have a more complete comparison, particularly for a quasi-equi-biaxial stress state. Other failure criteria could also be considered, for example, taking into account the shear stresses in the ceramic. Future developments of the software forecast the introduction of this type of criteria.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Fracture strength distribution of the initial and the damaged material.

Fig. 3 .

 3 Fig. 3. Weibull's plot of the four-point bending strength at 1423 K of ceramic specimens oxidised previously at 1423 K for 50 h.

Fig. 4 .

 4 Fig. 4. Stress±strain curves obtained from torsion test at 1423 K and calculated from the model. Fig. 5. Ceramic specimens and the workpiece used for the fracture tests.

Fig. 6 .

 6 Fig. 6. Orthoradial stress distribution obtained from the forging test simulation after a punch displacement of 4 mm.

Fig. 7 .

 7 Fig. 7. Experimental and calculated failure probability of the ceramic specimen.