Maxime Ployard
email: maxime.ployard@ec-lille.fr

Stéphane Brisset
email: stephane.brisset@ec-lille.fr

Florent Delhaye
email: florent.delhaye@ec-lille.fr

Julien Pouget
email: julien.pouget@sncf.fr

COMPARISON OF OPTIMIZATION FORMULATIONS TO DESIGN AN HYBRID RAILWAY POWER SUBSTATION

Keywords: Hybrid railway power substation, multi-sources system, energy storage, design, optimization, linear programming

This paper presents the design of a Hybrid Railway Power Substation by optimization. Two formulations are investigated to obtain the structural characteristics and the energy management strategy that minimize the total cost of HRPS during an operating period of one year. Energy formulation associated with sparse matrices allows optimizing on a wide time period with a short time step.

INTRODUCTION

Because of the rail network development, power demand of railway power substation (RPS) becomes increasingly important [START_REF] Aeberhard | Railway traction power supply from the state of the art to future trends[END_REF]. The integration of renewables energies could compensate this need and reduce pollution. Nevertheless the known intermittency factor of these energies makes difficult to estimate the amount of produced power [START_REF] Clastres | Ancillary services and optimal household energy management with photovoltaic production[END_REF]. By introducing a storage device, the hybrid RPS (HRPS) may operate without the grid and could ensure a good power quality [START_REF] Aeberhard | Railway traction power supply from the state of the art to future trends[END_REF]. The storage device should be accurately used according to the consumed power, the produced power and the evolution of prices fixed by power supplier [START_REF] Rigo-Mariani | A fast optimization strategy for power dispatching in a microgrid with storage[END_REF]. To enhance the design of HRPS, the energy management strategy should be suitable in the long term. However, a high number of unknowns is induced by these temporal variables and thereby two different formulations are investigated within the linear programing. By taking into account economical aspects, the optimal design of each source and the energy management strategy are known for a time step of ten minutes during one year.

1-MODEL AND OPTIMIZATION PROBLEM

The studied HRPS is presented in Fig. 1, it is composed of photovoltaic panels, wind turbines and a storage system connected to the substation through an AC bus. Renewable energy production are illustrated by one-way power flows on Fig. 1, depends of the photovoltaic panels surface 𝑆 𝑝𝑣 and the swept surface of wind turbines 𝑆 𝑤 . For the storage device, the quantity of maximal energy 𝐸 𝑚𝑎𝑥 defines the capability of storage, then power 𝑃 𝑠𝑡𝑜 (𝑡) or energy 𝐸 𝑠𝑡𝑜 (𝑡) can both express the temporal fluctuation of storage device for the optimization process. According to measurements of wind speed and solar irradiation and for a given consumption during one year, all power flows of this multi-source system have to satisfy the power balance given in (1) [START_REF] Pankovits | Design and operation optimization of a hybrid railway power substation[END_REF]. The power from the grid 𝑃 𝑔𝑟𝑖𝑑 is decomposed in two parts, the power received 𝑃 𝑔𝑟𝑖𝑑 + and the power sent 𝑃 𝑔𝑟𝑖𝑑 -.

(1)

𝑃 𝑔𝑟𝑖𝑑 (𝑡) = 𝑃 𝑔𝑟𝑖𝑑 + (𝑡) + 𝑃 𝑔𝑟𝑖𝑑
The objective function is expressed in (2) for an annual operating mode by considering investments cost of each device and also the purchasing 𝐶 𝑔𝑟𝑖𝑑 + and selling 𝐶 𝑔𝑟𝑖𝑑 of electricity [START_REF] Rigo-Mariani | A fast optimization strategy for power dispatching in a microgrid with storage[END_REF].

(2)

𝐶 = 𝛥𝑡. ∑ [𝐶 𝑔𝑟𝑖𝑑 + (𝑡) • 𝑃 𝑔𝑟𝑖𝑑 + (𝑡) + 𝐶 𝑔𝑟𝑖𝑑 -(𝑡) • 𝑃 𝑔𝑟𝑖𝑑 -(𝑡)] 𝑇 1
+ 𝐶 𝑠𝑡𝑜 𝐸 . 𝐸 𝑚𝑎𝑥 + 𝐶 𝑠𝑡𝑜 𝑃 . 𝑃 𝑚𝑎𝑥 + 𝐶 𝑃𝑉 . 𝑆 𝑃𝑉 + 𝐶 𝑊 . 𝑆 𝑊

In Table 1, two optimization problems are proposed through energy and power formulations of storage device. In power formulation, the stored energy is computed from the initial state of charge and also from the power stored at previous time steps. On the contrary, energy formulation considers only the maximal power of the storage. Dimensional variables of HRPS and grid power flow are the same on both optimizations.

2-COMPARAISON OF FORMULATIONS

For both optimization problems, linear programming is used. However it requires specific matrices to express equality and inequality constraints. Due to a time step of ten minutes, each temporal variable induces 52560 unknowns for one year cycle, therefore sparse matrices are employed. On Table 2, both formulations are compared according to the size of equality and inequality matrices and the memory required to store them. Fig. 2.a. and Fig. 2.b. present respectively the nonzero elements of equality and inequality matrices for power and energy formulations during one day. On Fig. 2.a., the power formulation induces a higher number of nonzero elements which makes difficult the storage of matrices for a high number of unknowns. By employing energy formulation, the memory required is greatly reduced and the optimization of HRPS is achievable in about one minute.

Figure 1 .

 1 Figure 1. Overview of Hybrid Railway Power Substation.

Figure 2 .

 2 Figure 2. Comparison of matrices for one day

Table 1 . Description of each optimization problem formulation Formulation power energy

 1

	Variables	common specific	𝑆 𝑃𝑉 , 𝑆 𝑊 , 𝐸 𝑚𝑎𝑥 , 𝑃 𝑚𝑎𝑥 , 𝐸 𝑖𝑛𝑖𝑡 , 𝑃 𝑔𝑟𝑖𝑑 + (𝑡), 𝑃 𝑔𝑟𝑖𝑑 -(𝑡) 𝑃 𝑠𝑡𝑜 (𝑡) 𝐸 𝑠𝑡𝑜 (𝑡)
		energy bounds	0 ≤ 𝐸 𝑖𝑛𝑖𝑡 -𝛥𝑡. ∑ 𝑃 𝑠𝑡𝑜 (𝑖) 𝑡 1	≤ 𝐸 𝑚𝑎𝑥	0 ≤ 𝐸 𝑠𝑡𝑜 (𝑡) ≤ 𝐸 𝑚𝑎𝑥
	Constraints	power bounds	-𝑃 𝑚𝑎𝑥 ≤ 𝑃 𝑠𝑡𝑜 (𝑡) ≤ 𝑃 𝑚𝑎𝑥	-𝑃 𝑚𝑎𝑥 . 𝛥𝑡 ≤ 𝐸 𝑠𝑡𝑜 (𝑡 -1) -𝐸 𝑠𝑡𝑜 (𝑡) ≤ 𝑃 𝑚𝑎𝑥 . 𝛥𝑡
		periodicity	∑ 𝑃 𝑠𝑡𝑜 (𝑖) = 0 𝑇 1		𝐸 𝑠𝑡𝑜 (𝑇) = 𝐸 𝑖𝑛𝑖𝑡

Table 2 . Comparison of memory requirements for both formulations

 2

	Formulation	Expression	power	Memory	Expression	energy	Memory
	Size of matrices	(5. 𝑇 + 1) x (3. 𝑇 + 5)	332 GB (full)	(4. 𝑇 + 2) x (3. 𝑇 + 5)	266 GB (full)
	Number of nonzero	14. 𝑇 + 𝑇 2		67 GB (sparse)	14. 𝑇 + 4		18 MB (sparse)


Raffaele D'elia
email: raffaele.delia176@gmail.com

Gérard Bernhart

Jouke Hijlkema

Thierry Cutard

Experimental analysis of SiC-based refractory concrete in hybrid rocket nozzles

Keywords:

Hybrid propulsion represents a good alternative to the more widely used liquid and solid systems. This technology combines some important specifications of the latters, as the possibility of re-ignition, thrust modulation, a higher specific impulse than solid systems, a greater simplicity and a lower cost than liquid systems. Nevertheless the highly oxidizing environment represents a major problem as regards the thermo-oxidation and ablative behavior of nozzle materials. The main goal of this research is to characterize a silicon carbide based micro-concrete with a maximum aggregates size of 800 μm, in a hybrid propulsion environment. The nozzle throat has to resist to a highly oxidizing polyethylene/nitrous oxide hybrid environment, under temperatures up to 2900 K. Three tests were performed on concrete-based nozzles in HERA Hybrid Rocket Motor (HRM) test bench at ONERA. Pressure chamber evolution and observations before and after tests are used to investigate the ablated surface at nozzle throat. Ablation behavior and crack generation are discussed and some improvements are proposed.

Introduction

Although hybrid propulsion has been introduced since many years, it is only in the last two decades that in view of cheap and reliable access to space has attracted the interest of many researchers and some companies. The environment produced by this kind of propulsion is more oxidizing than that of solid and liquid rockets, inducing a more pronounced ablation of nozzle's throat. This represents an important challenge as regard performances and nozzle's operation time. Usually graphite nozzles are used for hybrid rocket motors, without posing particular attention to this aspect. Authors observed a significant lack of works conducted in this field. Ablation was more studied in the case of solid rocket motors [START_REF] Maisonneuve | Ablation of solid-fuel booster nozzle materials[END_REF][START_REF] Kuo | Fundamental Understanding of Propellant/Nozzle Interaction for Rocket Nozzle Erosion Minimization Under Very High Pressure Conditions[END_REF][START_REF] Acharya | Effect of pressure and propellant composition on graphite rocket nozzle erosion rate[END_REF][START_REF] Bianchi | Modeling of Ablation Phenomena in Space Applications[END_REF][START_REF] Bianchi | Carbon-carbon nozzle erosion and shape change in fullscale solid-rocket motors[END_REF][START_REF] Bianchi | Chemical erosion of carbon-phenolic rocket nozzles with finite-rate surface chemistry[END_REF]. One of the few works about erosion in hybrid rocket nozzles is that of Nasuti and Bianchi [START_REF] Bianchi | Numerical analysis of nozzle material thermochemical erosion in hybrid rocket engines[END_REF]. They conducted a numerical study on the erosion of graphite nozzles, coupling Navier-Stokes equations with thermo-chemical ablative model, taking into account heterogeneous chemical reactions at nozzle surface, chemical diffusion speed in the boundary layer, interaction between ablation products and boundary layer and conduction in the material. In this frame, the French Space Agency (CNES) launched in 2005 the PERSEUS project, to develop new technologies applicable to the satellite launchers of the future. The project is mainly carried out by students in Master of Science and PhD courses. This work has started in 2009, in the context of the PERSEUS/PEGASE macro-project, leaded by AJSEP and engineers from Herakles Company. Its purpose is to design, develop and test a new nozzle, made of an SiC (Silicon Carbide) based micro-concrete. In 2011 a new SiC-based micro-concrete with a maximum aggregates size of 800 μm has been introduced, to face the problem of an important porosity degree in the material with a 3 mm aggregates size concrete. The use of concrete in nozzles is original, as in literature authors didn't find any work about this specific application. If the literature review is extended to concretes exposed to very high temperature and/or heat flux, some works can be found. It is the case of the Molten Core Concrete Interaction (MCCI), where liquid radioactive corium, generated by the fusion at around 2700 K of the nucleus of a nuclear reactor, interacts with refractory concrete that surround it [START_REF] Pakhomov | Fig. 11. HERA17 test: SEM images of longitudinal (a) and circumferential (b) breaking facets after test. It is possible to observe the concrete surface and the porous region. In figure (a)[END_REF][START_REF] Beshta | Interaction between a UO 2 þ x -ZrO 2 -Fe(Cr, Ni)O y melt and a concrete based on ZrO 2[END_REF][START_REF] Maruyama | A study on concrete degradation during molten core/concrete interactions[END_REF][START_REF] Journeau | European experiment on 2-D molten core concrete interaction : hecla and vulcano[END_REF]. Another case concerns the interaction between concretes and hot outlet gases in flame trench [START_REF] Liang | Final Report on the Study of Flame Trench Concrete Material[END_REF]. In a previous paper, D'Elia et al. [START_REF] D'elia | Preliminary tests of silicon carbide based concretes for hybrid rocket nozzles in a solar furnace[END_REF] studied the behavior of this concrete in standard atmosphere (0.84 bars of pressure), under concentrated solar flow up to 13.5 MW/m 2 , using a solar furnace. Surface erosion in these conditions is governed by silicon carbide active oxidation, with linear mass loss and a maximum erosion speed of around 5 μm/s. The goal of the present work is to extend this study to a hybrid propulsion environment, using the ONERA's HERA motor test bench. This study is mandatory because the environment in HRM is expected to be different with respect to the solar furnace conditions, in terms of pressure, flow velocity, oxidizing partial pressure and chemical composition. Erosion, ablation and oxidation phenomena observed in [START_REF] D'elia | Preliminary tests of silicon carbide based concretes for hybrid rocket nozzles in a solar furnace[END_REF] will be referred as Static Oxidation, in opposition to the Ablation phenomenon observed at nozzle throat, under an HRM environment.

Experiments

The hybrid rocket motor test bench used to perform fire tests on concrete based nozzles is named HERA (Fig. 1) and is located in Mauzac (FRANCE), at ONERA laboratory. It is composed of a polyethylene (PE) fuel cylinder, with a circular channel in the middle, held within an aluminum cylinder and mounted on a test bench table with steel supports. The channel within the PE fuel works as combustion chamber of the hybrid rocket motor. Nitrous oxide (N 2 O) is used as oxidizer. It is kept in liquid state in a 3.5 liters tank before use, under a 70-80 bars pressure (Fig. 1). Oxidizer is injected inside a pre-combustion chamber, where it is atomized, with a valve, located back to the combustion chamber. Oxidizer passes then through the combustion chamber: combustion is started with a pyrotechnic igniter, and can be stopped simply cutting off the oxidizer mass flow inside the chamber. A 3-6 bars chamber pressure is obtained for a 16 mm diameter of the nozzle's throat cross-section. These values are measured with two pressure sensors, located at the beginning and at the end of the combustion chamber. The second one is less reliable, because of its position close to the nozzle entrance. This is the reason why only data from the first sensor are considered in this paper.

Nozzles used in this study are made of a Silicon carbide based concrete (80.0 wt% SiC, 14.0 wt% Al 2 O 3 , 5.0 wt% SiO 2 , and 1.4 wt% CaO). They are designed using ONERA standard test geometry, with a contraction ratio (A c /A t ) of 10.5 and an expansion ratio (A e / A t ) of 4.0. ONERA CAD files were used to design mold geometry on CATIA. A resin-based mold was chosen to cast the concrete. Concrete's processing route is based on the same process explained in a previous work [START_REF] D'elia | Preliminary tests of silicon carbide based concretes for hybrid rocket nozzles in a solar furnace[END_REF]. Concrete is mixed with a 8.3 wt% water content during 6 min. Then it is poured into the mold, located on a vibrating table, during 15-20 min, to avoid macro-porosities generation. After that, material is maintained at 20 °C under plastic films during 48 h to ensure correct hydration. Concrete nozzle is then unmolded and stabilized at 110 °C during 24 h, to ensure water evaporation and to avoid material spalling at high temperature. Nozzle's samples are finally fired at 1200 °C during 5 h, to stabilize the microstructure, increase Young's modulus and reduce liquid phases generation during test at high temperature [START_REF] D'elia | Preliminary tests of silicon carbide based concretes for hybrid rocket nozzles in a solar furnace[END_REF][START_REF] Hewlett | Lea's Chemistry of Cement and Concrete[END_REF][START_REF] Lee | Castable refractory concretes[END_REF]. Material's properties are summarized in Table 1. After manufacturing (Fig. 2), nozzle is mounted at the end of the combustion chamber, directly in contact with the aluminum cylinder.

Three nozzles were tested in HERA motor, using ONERA geometry. The first two nozzles (HERA17 and 18) were only made with micro-concrete. The last one (HERA20) was made with micro-concrete, reinforced with 1 wt% of stainless steel short fibers (0.55 mm  25 mm) [START_REF] Cailleux | Study of a ceramic refractory reinforced with metallic fibres: from the microstructure to the mechanical behaviour[END_REF]. Processing route of the latter is slightly different from the classical method presented before. Fibers are first mixed with concrete during 5 min, before the addition of water. After this step, material's processing route is the same as for the others two nozzles.

Results

First nozzle was tested during a fire test labeled HERA17, with a combustion time of 20 s. The other two nozzles were tested during fire tests HERA18 and HERA20, with a combustion time of 10 s. Pressure curves are presented in Fig. 3 for the three nozzle tests. During HERA17, chamber pressure (Fig. 3b) increases quickly from 0.1 (atmospheric pressure) to 0.36 MPa in the first 0.7 test seconds. This phenomenon is generated by the transitory phase of ignition of the hybrid rocket motor. Then combustion starts and pressure increases from 0.36 MPa to 0.51 MPa after 9.2 s. After a stationary phase of 1.7 s, a pressure decrease is observed. This decrease continues until the end of combustion, after 20.2 s, with a pressure of 0.32 MPa. A second transitory phase of 0.5 s is observed, due to the motor shutdown: pressure decreases down to 0.06 MPa after 20.8 s, with a de-pressurization with regard to the atmospheric pressure, to reach again an equilibrium state after 25 s. This last phase of pressure equilibration is not considered in overall combustion time. Combustion of HERA17 motor can then be divided into 3 zones: ignition transitory, fuel regression and shutdown transitory. The end of ignition is defined as the moment when the injection pressure becomes constant (Fig. 3a). The end of combustion is the moment when the injection pressure falls very quickly, because no more oxidizer is injected in the combustion Fig. 1. HERA motor before a test. It is possible to notice the tank, the valve, the combustion chamber and the concrete nozzle. The igniter and the pre-chamber are not visible, because they are inside the engine.

Table 1

Physical and mechanical properties of the SiC-based micro-concrete [START_REF] D'elia | Preliminary tests of silicon carbide based concretes for hybrid rocket nozzles in a solar furnace[END_REF].

Property

Micro-concrete Density (kg/m 3 ) 2580 Coefficient of thermal expansion (10 À 6 K À 1 ) 5.75 Young's modulus (20-1200 °C, GPa) 75.0 Open porosity (vol%) 17.0 Fig. 2. SiC-based micro-concrete nozzle after manufacturing. chamber (Fig. 3a).

During fuel regression 3 sub-domains are observed: a. pressure increase, b. stationary pressure zone, c. pressure decrease. During this last phase, jet deviation and particles ejection are observed (Fig. 4a). Similar observations can be made for the test HERA18 and 20. In the last test HERA20, fuel regression seems to be almost completely stationary, after a little phase of pressure increase. According to the results of HERA17, these last two tests were performed with a combustion time of 10 s, to avoid jet deviation and particle ejection (Fig. 4b), and for a better understanding of ablation phenomenon. Durations of different phases of combustion in the three tests performed are summarized in Table 2.

To correctly evaluate material behavior, some additional information is necessary, as the oxidizing species molar fraction and the static temperature at nozzle throat. Open source Rocket Propulsion Analysis (RPA) software a is used to compute equilibrium calculation of the gas-dynamic flow in the nozzle. Combustion parameters are calculated introducing the following information into the software: chamber pressure, oxidizer and fuel chemical formula, mean mixing ratio, contraction and expansion ratio.

Chamber pressure used in the calculation is the value measured in the stationary pressure zone. Mixing ratio was calculated as follow:

∫ = ̇( ) O F m dt M / 1 t t o f i f
where O/F is the mean mixing ratio, ṁo is the oxidizer mass flow [kg/ s], measured at the inlet valve, M f is the fuel mass loss [kg], obtained from measurements before and after the test, t i and t f are the start time and the end time of combustion. Using these parameters, the information presented in Table 3 is obtained with RPA software. This program allows also the calculation of the gasdynamic field in a case of ideal combustion (stoichiometric mixing ratio, O/F sto ). The value obtained in this case is 7.46, for a chamber pressure ranging from 0.4 to 0.6 MPa.

Images of nozzle's throat cross-section are taken before and after test, with a telecentric lens (TCZR072, OPTO Engineering), mounted on a 5 Megapixel Black/White camera (AVT Pike F505B, ALLIED Vision Technologies). This system allows a correct view of nozzle throat, avoiding chromatic aberration, normally generated using classical lens that introduces some errors in cross-section measurement during image post-treatment. In order to measure total ablated surface, images are treated with ImageJ open-source software. Ablated surface is calculated as the difference between the two surfaces, measured before (S initial ) and after (S final ) the test; an equivalent ablated radius (r abl ) is calculated as:

π π = - = - ( ) r r r S S 2 abl final initial final initial
Cross-section evolutions for three different tests conditions are presented in Figs. 5b-d, in comparison with the untested nozzle (Fig. 5a). Ablated surface and equivalent ablated radius are presented in Table 4.

Discussion

Separation of gas dynamic field

Ablation was only observed for the first two nozzles tested at ONERA (HERA17 and 18). The last tested nozzle (HERA20) didn't show any erosion at nozzle throat. Separation of gas-dynamic flow could be a direct consequence of the ablation during the first test (HERA17), as observed in Fig. 4a. This test had a duration of 20 s and a jet flow deviation was observed after 10 s. The asymmetric ablation of the nozzle throat, clearly observed in Fig. 5b, changes the gas dynamic field, introducing an instability, which leads to the separation. According to experimental observations [START_REF] Sutton | Rocket Propulsion Elements[END_REF][START_REF] Nguyen | Nozzle separation prediction techniques and controlling techniques[END_REF], this phenomenon is normally generated in an over-expanded regime of the nozzle, when atmospheric pressure exceeds nozzle exit pressure by a factor of 2.5-4. Under these conditions, gas flow separation takes place, with recirculation zones in the divergent section and with instabilities. In this case, separation seems to be strictly related to nozzle's throat asymmetric ablation: when ablation is not observed at nozzle throat (HERA20), gas flow separation is not observed during the test. In the case of HERA18 nozzle, ablation was observed, but not jet deviation. It could be due to the fact that in this case the asymmetry of ablation at nozzle throat is not pronounced as in the case of HERA17. According to these observations it is possible to affirm that the separation phenomenon is a direct consequence of the asymmetric ablation of nozzle throat.

Parameters influencing throat ablation rate

Throat ablation rate is a strong function of mixture ratio and chamber pressure [START_REF] Bianchi | Numerical analysis of nozzle material thermochemical erosion in hybrid rocket engines[END_REF][START_REF] Sutton | Rocket Propulsion Elements[END_REF][START_REF] Bunker | Hybrid rocket motor nozzle material predictions and results[END_REF].

A change of the mixture ratio of the test typically affects the flame temperature (and the static temperature nozzle throat) and the molar mass of oxidizing species at throat. This effect is clearly visible considering the difference between the second and the third test (HERA18 and 20), which have the same duration (10 s), but an important discrepancy in terms of average mixture ratio (7.64 vs. 5.69): this variation strongly affects the static temperature of the two tests (respectively 2900 K and 2673 K) and molar mass of oxidizing species (respectively 17% mol and 15.6% mol for H 2 O; 9.3% mol and 5% mol for CO 2 ; 0.9% mol and 0.0% mol for O 2 ). These values explain very well the strong difference observed in terms of ablated radius: the nozzle tested during the HERA18 test bench is subjected to an heavier thermo-chemical load, which activate earlier the ablation mechanisms.

The second parameter that has an influence on ablation speed is the chamber pressure. This is one of the most important parameters that can affect the throat erosion rate [START_REF] Acharya | Effect of pressure and propellant composition on graphite rocket nozzle erosion rate[END_REF][START_REF] Bianchi | Carbon-carbon nozzle erosion and shape change in fullscale solid-rocket motors[END_REF][START_REF] Bianchi | Numerical analysis of nozzle material thermochemical erosion in hybrid rocket engines[END_REF][START_REF] Bianchi | Thermochemical erosion analysis for graphite/carbon-carbon rocket nozzles[END_REF][START_REF] Thakre | Chemical erosion of carbon-carbon/graphite nozzles in solid-propellant rocket motors[END_REF]. It is interesting to notice that in the case of HRMs the ablation rate increases linearly with chamber pressure [START_REF] Bianchi | Numerical analysis of nozzle material thermochemical erosion in hybrid rocket engines[END_REF], which is the same behavior already observed for solid rocket motors [START_REF] Acharya | Effect of pressure and propellant composition on graphite rocket nozzle erosion rate[END_REF][START_REF] Bianchi | Carbon-carbon nozzle erosion and shape change in fullscale solid-rocket motors[END_REF][START_REF] Bianchi | Thermochemical erosion analysis for graphite/carbon-carbon rocket nozzles[END_REF][START_REF] Thakre | Chemical erosion of carbon-carbon/graphite nozzles in solid-propellant rocket motors[END_REF]. This effect can partially explain the difference observed in the first two tests performed (HERA17 and 18), where a difference of 0.8 bars is observed in terms of chamber pressure (5.1 bars and 5.9 bars respectively) during the stationary phase. However, this difference of pressure is too small to explain alone the reduction of 60% of the ablated radius, observed for the second test HERA18. In this first two tests, the mixture ratio are very close and so the molar mass of the oxidizing species and the static temperatures (Table 3). The only parameter that can explain the discrepancy of the ablated radius is the ablation time (20 s for HERA17 and 10 s for HERA18), which plays a fundamental part in the degradation of the nozzle throat and in the activation of the separation phenomenon, as described above.

Chamber pressure vs. nozzle throat ablation

In order to analyze the evolution of the pressure curves presented in Fig. 3b, it is useful to introduce the mass flow rate equation for an ideal compressible gas: This equation shows that there are three main causes of the change in chamber pressure observed in Fig. 3b: mass flow rate ( ṁtot ), characteristic exhaust velocity (c*) and the nozzle throat cross-sectional area (A t ). In order to evaluate how these parameters changes with time and their impact on pressure curves, it is useful to analyze them according to HRM theory, as presented in The fuel surface regression rate, ṙ , where r is the radius of the port, is usually determined by the oxidizer mass flux, generally following the expression developed by Marxman and Gilbert [START_REF] Marxman | Turbulent boundary layer combustion in the hybrid rocket[END_REF] and Marxman [START_REF] Marxman | Boundary layer combustion in propulsion[END_REF] in the early 1960s:

̇= * → = * ̇( ) m p A c p c m A 3
̇= ̇= ( ) ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ r a m A a G 4 o x n n 0 0 0
where a 0 and n are empirical constants, ṁo is the oxidizer mass ] and G 0 is the oxidizer mass flux [kg/m 2 s]. The exponent n is typically between 0.5 and 0.8 [START_REF] Cantwell | Recent advances in hybrid propulsion[END_REF]. The fuel mass flow ( ṁf , [kg/s]) can be express as follow:

Particle es ejection

ρ π ̇= ⋅ ⋅ ⋅ ⋅ = ⋅ ⋅ ̇⋅ ( ) - m r D L K Lm D 5 f f x o n x n 1 2 where ρ π = ⋅ ⋅ ⋅ - K a 4 n f n 0
1 . The characteristic exhaust velocity can be expressed as follow:

* = Γ ( ) c RT M / 6 C
where R is the gas constant [8.3144 J/K mol], T C is the combustion temperature of the mixture [K], which strongly depends upon the mixture ratio O/F, M is the molar mass of the mixture [g/mol] and Γ is an adimensional parameter depending only on the adiabatic index γ:

γ γ Γ= + ( ) γ γ + ( -) ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ 2 1 7 1 2 1
Measurements, made on HERA motor test bench, have shown that oxidizer mass flow, after first ignition phase, does not change during test. According to this experimental observation and to Eqs. (3-7), and recalling that the exponent is generally in the range 0.5 on o0.8, the following considerations can be pointed out:

1. As the port diameter D x increases with time and the oxidizer mass flow ( ṁo ) is constant, the oxidizer mass flux (G 0 ) decreases with time; 2. The first consequence of point 1 is the decrease with time of the fuel surface regression rate ṙ (Eq. ( 4)); 3. The second consequence of point 1 is the decrease of the fuel mass flow (Eq. ( 5)). Note that for an exponent 0.5 on o0.8, the overall exponent of the port diameter D x is negative (1-2n); 4. The decrease of the fuel mass flow increases the mixing ratio O/ F for a constant oxidizer mass flow, with a consequent evolution of the characteristic velocity c * . This evolution depends upon the initial value of the O/F ratio. Typically, the flame temperature shows a maximum for the stoichiometric condition (O/F sto ) and tends to equally decrease for fuel-or oxidizer-rich mixtures [START_REF] Bianchi | Numerical analysis of nozzle material thermochemical erosion in hybrid rocket engines[END_REF]. The effect of this trend of the flame temperature is a maximum for the characteristic velocity at the mixture ratio (Eq. ( 6)); According to Eq. ( 3), if no ablation is observed (HERA20), the nozzle throat cross-sectional surface is constant during the test and there are only two components influencing the chamber pressure: the overall mass flow (which decreases with the combustion time) and the characteristic velocity. Since a constant pressure is observed during the HERA20 test, it is very likely that the characteristic velocity increases during the test, balancing the effect of the overall mass flow. This is possible if the O/F initial value is below the optimal value: it is exactly the case of HERA20, where an average value of 5.67 is obtained (Table 3), clearly below the stoichiometric value of 7.46.

When ablation is observed (HERA17 and 18), the nozzle throat cross-section A t is not constant anymore and its evolution will influence the chamber pressure value. A t can only increase with time, because of nozzle throat erosion, with a consequent decrease of chamber pressure. However, a decrease of chamber pressure can also be generated by the decrease of the characteristic velocity, if the mixing ratio O/F exceeds the stoichiometric value during the test, and by the decrease of the overall mass flow, as discussed above. In Fig. 3b and in Tables 2, 3 sub-domains are observed during the fuel regression phase, characterized by three different evolutions of chamber pressure: a. pressure increase zone, b. stationary pressure zone and c. pressure decrease zone: a. In the case of HERA17 and 18, the average mixing ratio (8.55 and 7.64 respectively) are above the stoichiometric value (7.46), but, since an increase in chamber pressure is observed in the initial phase of combustion, it is very likely that the O/F in these phases are below the optimal value, generating an increase of c * and a consequent increase of pressure. It is not possible to affirm if in this phase nozzle ablation takes place, because the increase of pressure generated by the characteristic velocity could very well overwhelm both the effects of chamber pressure decrease generated by ablation (if it is already active) and by the overall mass flow ( ̇) m tot ; b. In the second phase of combustion a stationary zone is observed. During this phase the decrease in chamber pressure generated by ṁtot and (possibly) by the ablation of nozzle throat is balanced by the characteristic velocity increase; c. In the third phase a decrease in chamber pressure is observed.

During this phase, the characteristic velocity has probably passed the stoichiometric value and it cannot balance anymore the others two effects generated by the ablation and by ṁtot , with a consequent decrease of the chamber pressure.

Ablation speed at nozzle throat is not constant, but depends on several parameters, as gas flow temperature and speed, oxidizing species molar fraction, surface roughness. It is not possible to correctly describe the evolution of all these parameters during the test. For this reason it is very difficult to understand when ablation starts exactly without any other data than the chamber pressure curves. However during the third sub-domain of the fuel regression phase of HERA17 test, a linear decrease of chamber pressure (clearly visible in Fig. suggests that nozzle ablation is already active: if the ablation phenomenon was activated during this phase, a further decrease in chamber pressure would change the slope of the curve. The separation of the gas dynamic field, observed for HERA17 after 10 s of tests, could be another element confirming the presence of material ablation in this phase. On the other hand, in the case of HERA18 the linear zone is very limited and the lack of separation does not help in the interpretation of results.

Material surface analysis

SEM images provide a better understanding of the material behavior at the microscopic level. Fig. 7 shows the nozzle after HERA17 test: the nozzle cross-section is presented on the left (a) and a SEM image of the nozzle throat is showed on the right (b). In the latter, nozzle throat can be divided in two zones: the first one (zone 1) on the top, where liquid phases are generated from micro-concrete at high temperature; the second one (zone 2), on the bottom of the image, shows a virgin surface, covered by some whiter particles.

A mixture of solidified liquid phases and silicon carbide aggregates composes the porous region on the top (zone 1). SiC aggregates at this temperature level (2900 K) undergo to active oxidation and sublimation. All these phenomena increase the erosion at material surface, with regard to static oxidation conditions in solar furnace [START_REF] D'elia | Preliminary tests of silicon carbide based concretes for hybrid rocket nozzles in a solar furnace[END_REF]. However it is not possible to confirm active oxidation and sublimation of SiC, because no direct observation of these phenomena was possible. Liquid phases, generated during the interaction between concrete surface and gas dynamic flow, can be seen as a moving front on the concrete, from the convergent to the divergent side, following the direction of the supersonic flow. Jet separation from the surface of the concrete may have occurred at the limit of this region (zone 1): this could explain why in zone 2, concrete does not seem to be attacked by the high reactive flow. Fig. 8 shows a detail of zone 1, observed on the top of Fig. 7b. In a previous work [START_REF] D'elia | Preliminary tests of silicon carbide based concretes for hybrid rocket nozzles in a solar furnace[END_REF], surface of micro-concrete appears as completely different after exposition to a 13.5 MW/m 2 solar flow, with silicon-carbide aggregates that undergo active oxidation. A silicon carbide aggregate, during active oxidation (Fig. 9a andb), shows a classical oxidation crown, which is the consequence of SiO production at high temperature (42100 K). Gas static temperature at nozzle throat is 2900 K, but oxidation crown neither silicon carbide aggregates can be directly observed. In contrary, the surface is covered by a solidified eutectic phase, generated at high temperature and solidified during nozzle cooling.

These phases have surrounded SiC aggregates that are not visible by SEM; the interaction with supersonic gas flow induces a movement towards the nozzle exit, with an increase of the nozzle throat cross-section. Therefore, it is not possible to confirm that

D x D t L PE PE PE ( ρ f )
Combustion chamber Nozzle Fig. 6. HRM motor scheme.

carbide active oxidation plays an active role in nozzle ablation. However when combustion time exceeds 10 s, it is possible to observe particles ejected from the nozzle. These particles, which have redeposited on the virgin concrete surface (Fig. 7b), are probably generated by an intense active oxidation of silicon carbide at high temperature or by SiC sublimation. The dynamic ablation is therefore mainly generated by the interaction between the high temperature supersonic flow and the eutectic liquid phases, produced from the cementitious matrix at high temperature. The role of silicon carbide active oxidation/sublimation cannot be confirmed in this work.

It is interesting to notice that jet deviation, observed in Fig. 4a, is directed upwards; this observation might suggest a role of gravity in the liquid phases mobility at nozzle throat: the combined action of supersonic flow and gravity could generate liquid phases motion directed towards the exit of the nozzle and downwards, deviating the flow upwards. To validate this hypothesis it could be interesting to perform a firing test with the nozzle axis perpendicular to the ground (vertical test), since in this configuration the gas flow and the gravity would have the same direction, avoiding the asymmetry of the system.

After test, the first two nozzles (HERA17 and 18) showed two different types of crack: a longitudinal one (Fig. 10a, type 1) and a circumferential one (Fig. 10a andb, type 2). The third nozzle (HERA20) did not show any crack (Fig. 10c). The first two nozzles are composed of micro-concrete without any reinforcement. The last one contains 1 wt% of short metallic fibers. HERA17 test duration was of 20 s and jet deviation was observed, while HERA18 test duration was of 10 s, without jet deviation. This difference appears to be critical in mechanical behavior of microconcrete, since a longitudinal crack appears when jet deviation is observed. The asymmetry of the gas dynamic field might suggest a tangential thermal gradient inside the structure, which could lead to this specific type of crack. HERA18 and 20 tests confirm the connection between jet deviation and the first kind of crack.

Otherwise, the second type of crack seems to be connected with the system configuration: the assembly/disassembly of the nozzle, its cooling phase or the geometry of the structure could induce some external stresses in the generating this kind of crack. The Fig. 10b might suggest that this crack was initiated during the combustion or during the assembly, because of the polyethylene black deposit on the crack surface. According to these observations it is possible to exclude that this second type of crack is initiated during the disassembly of the system or during material cooling, but it is not possible to confirm precisely if it is induced by the high thermal gradient or by stresses generated during the assembly phase. The two types of crack are the result of quasi-brittle rupture processes. Cementitious matrix boiling and/ or in-depth SiC aggregates oxidation/sublimation mechanisms at high temperature (T 42100 K) generate micro-and macro-porosities beneath the concrete surface. The result is a thin layer of porosities in depth in the material. The thickness of this layer increases with the diffusion of the oxidizing species, which depends on external conditions of temperature and pressure, but also on surface condition. The diffusion of these species in the material at high temperature contributes to material degradation and porosities generation. SEM images of the two breaking facets (vein profile side) are shown in Fig. It is possible to notice two regions in both of the crack type (1 and 2 in Fig. 10a): a first zone of higher porosity, which reaches a depth close to 700 μm for the longitudinal crack and close to 100 μm for the circumferential crack. This strong difference likely depends on the surface condition of the zone where the crack is observed: the longitudinal crack (Fig. 11a) is observed beneath a very porous region of the concrete (zone 1 of Fig. 7b), where species diffusion is faster; the circumferential crack (Fig. 11b) is observed under a very compact region (zone 2 of Fig. where diffusion is more difficult.

This region represents the weakest zone in concrete, where it is most likely that fracture is initiated. Unfortunately these hypotheses need more studies to be confirmed. Therefore according to present results, it is not possible to confirm precisely the reasons of crack initiation and propagation. However the introduction of short metallic fibers in the third nozzle has shown the possibility to avoid crack initiation in this material, if combustion time is lower than 10 s. This has been observed during the last test on concrete-based nozzle (HERA20), where no crack was observed after the test. Metallic fibers are randomly oriented inside the cementitious matrix because of vibration during concrete pouring into molds. As shown in previous works [START_REF] Cailleux | Study of a ceramic refractory reinforced with metallic fibres: from the microstructure to the mechanical behaviour[END_REF][START_REF] Cailleux | Pullout of metallic fibres from a ceramic refractory matrix[END_REF][START_REF] Cailleux | Pullout of steel fibres from a refractory castable: experiment and modelling[END_REF][START_REF] Orteu | Application of stereovision to the mechanical characterisation of ceramic refractories reinforced with metallic fibres[END_REF], the introduction of short metallic fibers in a refractory concrete allows to avoid its brittle behavior. On the one hand, diffuse damage develops in the reinforced refractory because of stress fields induced during the firing cycle, by a different dilatometric behavior between fibers and concrete constituents. Damage deals with micro- in the cementitious matrix and with fiber/concrete debonding. On the other hand, fiber/concrete mechanisms take place in the material when mechanical loadings are applied either as an external loading or as loading due to thermal gradient. They deal with micro-crack and macro-crack bridging, with fiber/concrete debonding and interface phenomena and with fiber pull-out mechanisms. These various mechanisms make the refractory concrete more tolerant to damage and increase the non-linearity of its mechanical behavior. As a consequence, the reinforced concrete is able to develop a better resistance to deformation even if its stress capability is reduced.

Conclusions

Three SiC-based refractory micro-concrete nozzles have been tested in a hybrid rocket motor, under a highly oxidizing polyethylene/nitrous oxide environment. Nozzles were tested during 20 s (HERA17 test) and 10 s (HERA18 and 20 tests), under a static pressure at nozzle throat of around 3 bars and a static temperature at nozzle throat of 2700 K (HERA20 test) and 2900 K (HERA17 and 18 tests). Three combustion domains can be observed on pressure curves: ignition transitory, fuel regression and shutdown transitory. Image analysis of nozzle throat cross-section before and after test, gave an estimation of the ablated surface. This phenomenon was detected only in the first two nozzles (HERA17 and 18 tests), where the higher combustion temperatures and the more important concentration of the oxidizing species have probably played an important role in the activation of nozzle ablation. Separation of gas dynamic flow has been observed for a combustion time greater than 10 s in the case of HERA17. This phenomenon is strictly related to an asymmetric ablation at nozzle throat, generated by liquid phases interaction with supersonic gas dynamic flow. Gravity has probably an important role with regard to the asymmetric ablation, but a vertical test would be necessary to confirm it. No separation was observed in the case of HERA18, where the combustion time was probably too short to generate an asymmetric ablation of the nozzle throat. Silicon carbide active oxidation and sublimation was not observed directly on the attacked surface after test, but particle ejection suggests its presence during this test phase. Authors observed the strong dependence of ablation phenomenon upon static temperature and molar mass of oxidizing species at nozzle's throat and upon combustion time: an increase of these parameters increases the ablated surface. Two types of crack were observed: a first longitudinal one (type 1) is generated by jet separation and its consequent tangential thermal gradient; a circumferential one (type 2) depends on the stresses, generated during combustion or motor assembly. More tests would be necessary to better understand cracks initiation processes, propagation and their dependence on experimental conditions. Short metallic fibers addition avoids crack propagation in concrete, improving the durability of the nozzle at high temperature.

According to this study, SiC-based micro-concrete nozzle, with short metallic fibers addition, has a good thermo-mechanical resistance under a hybrid rocket motor environment. This material, therefore, could be a good candidate as a cheaper alternative to the classical nozzle materials, as graphite and refractory metals, to use in hybrid rocket motor. A necessary condition to its use might be the application of a specific surface treatment before the test, reducing surface porosities and improving its resistance to ablation, oxidation/sublimation and liquid phases generation. A possible approach could be the Chemical Vapor Deposition technique, to create a thin layer of less reactive material; another could be the use of a thin layer of thermoset ablative material, with a higher ablation heat, to delay ablation initiation.

Fig. 3 .

 3 Fig. 3. Injection (a) and chamber (b) pressures evolution versus time, during nozzle tests HERA17, 18 and 20 (chamber pressure measured at the beginning of the combustion chamber).

  represents the chamber pressure [Pa], * c is the char- acteristic exhaust velocity [m/s], ̇= ̇+ ṁ m m tot o f is the overall mass flow [kg/s] and A t is the nozzle throat cross-sectional surface [m 2 ].

Fig. 6 .

 6 In this figure, L is PE block/combustion chamber length [m], D x is the combustion chamber diameter [m], D t is the nozzle throat diameter [m], ṙ is fuel regression rate [m/s] and ρ f is fuel density [kg/m 3 ].

Fig. 4 .

 4 Fig. 4. (a) HERA17: jet deviation and particle ejection (white arrows) for a combustion time greater than 10 s. (b) HERA18: for combustion time lower than 10 s neither jet deviation nor particle ejection are observed. The same behavior is observed for HERA20.

Fig. 5 .

 5 Fig. 5. Nozzle throat cross-sections before (a) and after (b-d) test. (a) Untested nozzle. (b) HERA17: 20 s. (c) HERA18: 10 s. (d) HERA20: 10 s.

5 .

 5 The decrease of the fuel mass flow decreases the overall mass flow ( ̇= ̇+ ṁ m m tot f o ) for a constant oxidizer mass flow.

Fig. 7 .

 7 Fig. 7. HERA17 nozzle after test: (a) nozzle cross-section and (b) SEM image of the nozzle throat surface.

Fig. 8 .

 8 Fig. 8. HERA17 nozzle after test: SEM image of nozzle throat surface in the convergent side (zone 1).

Fig. 9 .

 9 Fig. 9. Micro-concrete tested in a solar furnace [13], under a concentrated solar flux of 13.5 MW/m 2 and a calculated surface temperature of 2400 K: (a) top view and (b) cross-section of SiC aggregates under active oxidation.

Fig. 10 .

 10 Fig. 10. HERA17 (a), HERA18 (b) and HERA20 (c) nozzles after test. Two types of cracks observed after HERA17 and HERA18 tests: longitudinal cracks (type 1) and circumferential crack (type 2).

Table 2

 2 Durations (Δt) and pressure drop (Δp) of the different combustion phases for the three tests performed at ONERA.

		HERA17		HERA18		HERA20
	Ignition transitory Fuel regression (p increase)	þ0.26 þ0.15	0.7 8.5	þ0.35 þ0.14	0.7 5.0	þ0.42 þ0.05	0.8 1.3
	Fuel regression (sta-	-	1.7	-	0.9	-	8.0
	tionary p)
	Fuel regression (p decrease)	À 0.20	9.3	À 0.06	3.4	-	-
	Shutdown transitory Overall combustion	À 0.26	0.5 20.7	À 0.49	0.9 10.9	À 0.48	0.5 10.6
	time

Δp (MPa) Δt (s) Δp (MPa) Δt (s) Δp (MPa) Δt (s)

Table 3

 3 Combustion parameters of the three tests performed at ONERA (* data are calculated with RPA software).

		HERA17	HERA18	HERA20
	p c (MPa)	0.51	0.59	0.57
	Average ṁo (g/s)	74.7	78.4	65.3
	M f (g)	175	103	109
	O/F	8.55	7.64	5.69
	p t (MPa) *	0.29	0.32	0.32
	p exit (MPa) * c * (m/s) *	0.03 1536	0.03 1552	0.02 1543
	M exit *	2.51	2.53	2.67
	T t (K) *	2897	2902	2673
	H 2 O (%mol) *	16.3	17.0	15.6
	CO 2 (%mol) *	10.3	9.3	5.0
	O 2 (%mol) *	2.0	0.9	0.0

Table 4

 4 Ablated surfaces and equivalent ablated radius of HRM tests.

	Test	Ablated surface (mm 2 )	Equivalent ablated radius (mm)
	HERA17	35.43	0.67
	HERA18	12.95	0.25
	HERA20	0.00	0.00

Acknowledgments

Authors would like to thank CNES and Midi-Pyrenees region for the financial support to the Ph.D. thesis. Special thanks to PERSEUS-PEGASE project leaders Jean Oswald and Christophe Magnière.


Gilles Dusserre
email: dusserre@enstimac.fr

Gilles Dour

Gérard Bernhart

In-situ determination of the heat flux density at the glass/mould interface during a glass pressing production cycle

Keywords: Glass pressing process, Heat flux density, Temperature measurement, Inverse method, FEM

The glass pressing process involves heat transfer between the glass gob and the forming tool which are among the most important parameters influencing the thermo-mechanical stresses in the moulds. The present paper presents the development of the instrumentation of a mould for the measurement of temperatures during the production cycle. These measurements are exploited with an inverse method to evaluate the heat flux densities at the working surface of the mould. The influence of each process stage and of the location at the surface of the mould on the thermal loadings are described. The evaluated heat flux densities are used as boundary conditions in a finite element calculation. The validity of these results are discussed taking into account the differences between experiment and calculation, the hypothesis of the inverse method and the time response of the thermocouples.

Introduction

Problematic and objectives

Glass pressing is a forming process in which a glass gob is squeezed between two tools, currently made out of steel, to produce hollow glass parts for chemical or cooking applications (see Fig. 1). The principal functions of the moulds are to give their shape to the glass and, in particular for the lower one, to extract heat from the glass. This heat exchange between gob and mould is all the more important that the dilatation gradient created by this way generates stresses in the mould. It has been demonstrated in a previous paper [START_REF] Dusserre | Thermo-mechanical stresses in cast steel dies during glass pressing process[END_REF] that thermal stresses are the most important stresses, compared to mechanical stresses such as those resulting from the glass compression. In regard to the limited life time of the moulds made of cast tempered martensitic stainless steel and the cost of such tools [START_REF] Pöting | Future trends in hollow glass production[END_REF], it is very important for the glass industry to quantify the heat flux density exchanged between the gob and the mould in order to optimize the design, the material and the use of the moulds. It is also important to optimize the productivity of the process.

The final aim of our industrial research project being to calculate the thermo-mechanical stresses in the mould and considering the difficulties to measure the surface temperature of glass because of its semi-transparent properties [START_REF] Linhart | Temperature measurement in the glass melt[END_REF], it has been decided to focus our attention on the heat flux density, which could be used as a boundary condition in a finite element simulation, instead of the heat transfer coefficient. To understand how the process conditions affect the heat exchange between the gob and the mould, it is here proposed to measure temperatures in the lower mould in use during the industrial production cycle. The aim of this work is to evaluate the heat flux density exchanged at the surface of the die, using an inverse method to exploit the temperature measurements.

The paper proposes to explain our instrumentation strategy, the results in terms of temperature and heat flux density as a function of the die location and some process parameters. Finally the precision of the results will be discussed with regards to the sensors and the inverse method used.

Nomenclature

Bibliographic background 1.2.1. Heat exchange between two bodies

The heat transfer at the interface between two materials is controlled by material physical and mechanical properties [START_REF] Rawson | Comments on "Heat transfer modelling of the parison forming in glass manufacturing[END_REF], the surface roughness [START_REF] Höhne | Measurement and mathematical modeling of the heat transfer in the glass forming process, in consideration of the heat transfer coefficients and radiation influences[END_REF] and the contact pressure [START_REF] Storck | Heat transfer modelling of the parison forming in glass manufacturing[END_REF]. Due to the roughness of the surfaces, the true contact zone is not a continuum but is constituted by local contacts and local gaps alternatively [START_REF] Loulou | Estimation de la résistance thermique de contact durant la solidification du verre[END_REF]. The temperature of the two surfaces in contact is delicate to measure if not to define. The usual practice is to define a temperature at a distance from the surface profile, i.e. ten times the height of the surface. From that definition of the surface temperature results that the two temperature for the two materials in contact are necessarily different [START_REF] Khachfe | Estimation de la résistance thermique de contact entre un métal et une pastille de verre[END_REF]. This difference can be related to the micro-constriction which is the origin of a thermal contact resistance. The heat transfer coefficient is defined as the inverse of the thermal contact resistance.

Empirical models have been proposed [START_REF] Storck | Heat transfer modelling of the parison forming in glass manufacturing[END_REF][START_REF] Viskanta | Analysis of heat transfer during glass forming[END_REF], for example by Pchelyakov and Guloyan [START_REF] Pchelyakov | Heat transfer at the glass mould interface[END_REF] which considers that an isolating gap is present between the two opposite bodies. The heat transfer coefficient h is supposed to be the ratio between the conductivity of this gap k and its thickness δ (Eq. ( 1)).

h = k δ (1)
This model has been applied to glass pressing process describing the gap thickness at each stage of the process, taking into account the contraction of the glass, the contact pressure and some empirical parameters (see Eqs. ( 2)-( 4)).

-Before pressing:

δ = δ 0 + δ c (2)
-while pressing:

δ = δ P = C P + δ min (3)
-after pressing:

δ = δ end P + δ c ( 4
)
where δ c is calculated by integration from the thermal dilatation coefficient and the temperature field in the bulk of glass.

The limits of such model are the ignorance of true contact zones between the two opposite bodies. Then the value of the gap conductivity cannot be the air or lubricant one. The surface roughness of the mould and the surface tension of the glass are also ignored, whereas their value actually control the micro-constriction effect. The linear evolution of the heat exchange coefficient with contact pressure discussed in [START_REF] Keijts | Heat transfer in glass container production during the final blow[END_REF] for glass blowing is compatible with this model if the value of δ min could be neglected compared with the ratio between C and P . This condition is probably relevant for blowing process where contact pressure is not so high as in pressing process, but is not expected to apply to glass pressing.

Radiative heat transfer also may play an important role in the heat exchange between two bodies in contact, on one hand because each surface radiates the opposite one, and on the other hand, for glasses, because the whole bulk of a semi-transparent media radiates the surface of the opposite body (or the bulk if both bodies are semi-transparent media). The surface contribution is often taken into account by a simple expression using the emissivity of the surface, its temperature at a power 4 and the Stefan-Boltzmann constant [START_REF] Storck | Heat transfer modelling of the parison forming in glass manufacturing[END_REF][START_REF] Loulou | Estimation de la résistance thermique de contact durant la solidification du verre[END_REF]. The bulk contribution is more difficult to evaluate, for example solving the radiative transfer equation [START_REF] Khachfe | Estimation de la résistance thermique de contact entre un métal et une pastille de verre[END_REF][START_REF] Viskanta | Analysis of heat transfer during glass forming[END_REF][START_REF] Tarakanov | Some problems of precise measurements of heat transfer coefficients in glass melts. I: Measurements of effective conductivity[END_REF] and is rarely taken into account.

Heat transfer inverse problems

The lack of proper modeling of the micro-constriction thermal resistance and of the radiative heat transfer have pushed investigators to determine surface temperature and interfacial heat flow from bulk measurements. The difficulty here is the extrapolation to the surface. Historically speaking, Beck was the first to propose a rigorous method to extrapolate surface information from bulk measurements: inverse problem method.

Inverse methods have then been widened to determine some characteristics of heat transfer such as materials properties [START_REF] Sawaf | An inverse analysis to estimate linearly temperature dependent thermal conductivity components and heat capacity of an orthotropic medium[END_REF], heat sources [START_REF] Abou Khachfe | Résolution numérique d'un problème d'estimation de source thermodépendante dans un domaine bidimensionnel[END_REF] or heat transfer coefficient [START_REF] Abou Khachfe | Determination of heat sources and heat transfer coefficient for two-dimensional heat flow: Numerical and experimental study[END_REF]. The problem requires to combine a method to solve the direct problem, for example the thermal quadrupoles [START_REF] Maillet | Thermal Quadrupoles -Solving the Heat Equation through Integral Transforms[END_REF] or the finite differences [START_REF] Moreau | Méthodologie expérimentale pour la détermination du coefficient de transfert thermique à l'interface verre/outil[END_REF] in one dimension, the finite elements [START_REF] Abou Khachfe | Numerical solution of 2-d nonlinear heat conduction problems using finite-element techniques[END_REF] in two or three dimensions [START_REF] Jarny | A general optimization method using adjoint equation for solving multidimensional inverse heat conduction[END_REF], and an optimization technique, often iterative, for example the Beck's method [START_REF] Beck | Parameter Estimation in Engineering and Science[END_REF], the conjugate gradient [START_REF] Khachfe | Résolution numérique de problèmes inverses 2D non linéaires de conduction de la chaleur par la méthode des éléments finis et l'algorithme du gradient conjugué -Validation expérimentale[END_REF] or the Levenberg-Marquart method [START_REF] Moreau | Méthodologie expérimentale pour la détermination du coefficient de transfert thermique à l'interface verre/outil[END_REF]. The convergence of the optimization algorithms is difficult to ensure because of the high sensibility to the initial conditions, to the physical and numerical parameters of the direct model and to the choice of the optimization objectives. The advantage of the thermal quadrupoles method is that it does not require any spatial discretization, excluding any numerical diffusion phenomenon. It also allows a limitation of the time calculation because it is not necessary to evaluate the whole temperature field. However, its most important inconvenient is the impossibility to take the thermal variation of the material parameters into account, unlike discrete methods.

Experimental heat transfer determination

In literature, experiments have been conducted to determine the value of the heat transfer coefficient at the contact between glass and mould. They all consist in using some temperature measurements to evaluate by inverse method the value of the heat flux density. This value is divided by the difference between the two surface temperature to obtain the heat transfer coefficient. The results of McGraw et al. [START_REF] Mcgraw | Transfer of heat in glass during forming[END_REF][START_REF] Brauns | Finite element calculation of thermomechanical stresses in candidate ceramic components for press and blow moulds[END_REF] have been modeled by Jones [START_REF] Jones | BGIRA research report n • 29[END_REF] in terms of time evolution by a simple exponential expression. Moreau [START_REF] Moreau | Méthodologie expérimentale pour la détermination du coefficient de transfert thermique à l'interface verre/outil[END_REF] made measurements of temperature in a substrate in contact with glass to evaluate by inverse method (Levenberg-Marquart) the heat transfer coefficient. Loulou [START_REF] Loulou | Estimation de la résistance thermique de contact durant la solidification du verre[END_REF] did similar experiments to study the radiation effect on the thermal contact resistance using a one-dimensional modeling of the radiative transfer equation. Some authors proposed a constant value in conditions close to those of the glass pressing process (about 8 kW m -2 K -1 for Viscanta and Lim [START_REF] Viskanta | Analysis of heat transfer during glass forming[END_REF], about 3 kW m -2 K -1 for Höhne et al. [START_REF] Höhne | Measurement and mathematical modeling of the heat transfer in the glass forming process, in consideration of the heat transfer coefficients and radiation influences[END_REF]). In fact, most authors have shown that this value decreases when the time of contact increases, from a maximal value between 9 and 13 kW m -2 K -1 to a minimal value between 1 and 3 kW m -2 K -1 after a few seconds [START_REF] Pchelyakov | Heat transfer at the glass mould interface[END_REF][START_REF] Moreau | Méthodologie expérimentale pour la détermination du coefficient de transfert thermique à l'interface verre/outil[END_REF][START_REF] Mcgraw | Transfer of heat in glass during forming[END_REF]. The contraction of the glass while cooling seems to be the origin of such a phenomenon. Even if the maximal value is of the same order in most cases, the comparison of the different results of the literature shows the importance of the experimental conditions on the evaluated heat flux density. That is why it is here proposed to evaluate the heat exchange at the surface of the mould in the actual industrial context.

Experimental setup and inverse method

A technological know-how has been developed by CROMeP to investigate heat transfer occurring in forming processes using the Beck's method. In particular, the high pressure die casting process has been studied [START_REF] Dour | Development of a nonintrusive heat transfer coefficient gauge and its application to high pressure die casting: Effect of the process parameters[END_REF] from thermocouple and pyrometer measurements.

The present industrial process involves sixteen moulds positioned on a carousel and located successively at each stage of the process (see Fig. 2). The rotation cycle time is 48 s. The constraints for the data acquisition system resulting from such industrial context are:

-a high ambient temperature of about 60 • C, -a reduced size, -the impossibility of any physical link with the ground during the acquisition because of the rotation of the carousel, -the necessity of a quick assembly of the mould after its furnace preheating.

The moulds are made of GX30Cr13 cast stainless steel (see Table 1). Before machining and polishing, the raw state is obtained by sand casting followed by a heat treatment consisting an austenitization, a quenching and a tempering (see Table 2). This heat treatment leads to a tempered martensitic microstructure with a limited hardness compared with traditional mould materials.

Experimental setup

The gauge designed especially for this application is based on the model developed by Dour [START_REF] Dour | Development of a nonintrusive heat transfer coefficient gauge and its application to high pressure die casting: Effect of the process parameters[END_REF] for the pressure die casting process. It is constituted by a body in which 3 housings of 0.3 mm diameter have been drilled at 1, 7 and 13 mm from the surface of the gauge. In these housings are stuck the hot junctions of 3 type K thermocouples diameter 0.25 mm (their thinness allows for a very short response time) for the measurements of the temperatures required in the inverse method. The temperature measured at 13 mm is imposed, the one at 1 mm is used for the optimization and the one at 7 mm is compared to the calculation to validate the results. The bodies of the gauges are machined from a mould that had been cast especially for that use, in order to be made out of the same materials (composition and solidification structure) than any other one. The gauges are covered with a jacket to protect the thermocouples and placed in an industrial mould in which 3 housings were machined (see Fig. 3).

This so instrumented mould being used in an industrial context, an unattended acquisition device has been made up to answer to the industrial constraints. An embedded power supply and data storage were required. The inverse method imposes a minimal acquisition frequency and a good accuracy of the measurements. These specifications have been achieved with an industrial automaton National Instrument cFP 2000 realizing the data acquisition at 30 Hz during the preheating and a few tens of cycles. The cold junction compensation of the thermocouples has been taken into account by measuring the temperature in the box containing the acquisition material with a PT100 sensor. Each signal providing from a thermocouple has been amplified by an operational amplifier to obtain a resolution lower than 0.1 • C (low cost but time consuming calibration for each thermocouple).

Design of the thermal gauge

The inverse method that has been used in the present study is founded on the thermal quadrupoles formulation [START_REF] Maillet | Thermal Quadrupoles -Solving the Heat Equation through Integral Transforms[END_REF] (Laplace transformation of the heat conduction equation in one dimension) and the optimization problem is solved using the Beck's method [START_REF] Beck | Parameter Estimation in Engineering and Science[END_REF] as described by Broucaret [26].

This method requires two temperature measurements: the first one is imposed as a boundary condition in the direct problem and the second one is used to minimize the differences between experiment and calculation. A relevance diagram (Fig. 4) has been mapped out [START_REF] Dour | Development of a nonintrusive heat transfer coefficient gauge and its application to high pressure die casting: Effect of the process parameters[END_REF][START_REF] Dour | Recommendations and guidelines for the performance of accurate heat transfer measurements in rapid forming processes[END_REF] to determine the validity of the results as a function of some dimensionless parameters relevant to the acquisition (sampling period dt, distance between the surface and the first measurement location e), to the process (characteristic time of the heat exchange τ , chosen equal to 1 s a priori, corresponding to the duration of the application of the pressure on the gob) and to the mould material (diffusivity κ, 4.9 × 10 -6 m 2 s -1 for the studied steel). In order to be located in the middle of the validity domain, it has been retained a distance of 1 mm between the first measurement location and the surface, and a sampling time of 0.03 s.

Results and discussion

The instrumented mould presented above has been used in the factories of the group Arc International Cookware at Sunderland (UK) and Châteauroux (France). The results presented here concern the measurements realized at Sunderland. The preheating temperature is about 370 • C and about twenty cycles are required to stabilize the mean temperature. In a stabilized cycle, the surface temperature varies between 430 and 605 • C.

Accuracy of the heat flux estimation

From temperature measurement the heat flux and the surface temperature are estimated with the inverse method. The map of relevance (see Fig. 4) showed that the effectiveness of the inverse method is not straightforward. To map out this relevance diagram, the response of the direct problem to a heat flux step function has been used as data of the inverse method. The flux obtained for many values of the inverse method parameters have been compared to the initial heat flux step function to build the different zones of the diagram. The reference [START_REF] Dour | Recommendations and guidelines for the performance of accurate heat transfer measurements in rapid forming processes[END_REF] shows that using this method, no significant influence of the ntf parameter could be observed in the acceptable zone. We have chosen the parameters of the inverse method ( t) and of the sensors (e) so that it has the best chance to work according to our indications about the process (τ = 1 s).

Nevertheless with ntf = 1, the inverse method is perfectly unstable (the number of future instant, ntf, is the number of consecutive time step used by the Beck's method to evaluate the least squares error whose minimization allows the estimation of the heat flux density value at the next time step). In the acceptable zone, the results of inverse method should not depend on the ntf chosen. In the past experiences of Hamasaiid [START_REF] Hamasaiid | Etude des phénomènes thermiques à l'interface pièce/outil en fonderie[END_REF], the procedure indeed proposes to choose optimum ntf so that the noise in flux signal is as small as possible, but provided that maximum heat flux density remains in the range of maximum heat flux ± noise 2 with lower ntfs. Unluckily this procedure does not work for us. The effect of ntf is summarized in Fig. 5. We believe that the extremely brief duration of the peak heat flux (about 0.1 s instead of 1 s evaluated a priori due to the quickly decreasing heat flux curve) makes that the operative point on the map of relevance may be translated towards under or over evaluation or non-sense zone. Then it can be understood that the peak value depends on ntf. Nevertheless it has no influence on the maximum temperature. Considering that q ≈ k T

x ≈ 25 40 10 -3 = 10 6 W m -2 , a value as high as 6 MW m -2 does not seem reasonable.

An other way to evaluate the maximum value of the heat flux density is the normalized approach detailed in [START_REF] Dour | Thermal stresses and distortion in die casting processes: a new normalized approach[END_REF] if some information about the process are available, as the geometry, the material of the mould, the temperature variation during one , what is between the two values found with ntf equal to 3 and 4 (see Fig. 5). This method is very sensitive to the value of the heat flux duration and a little variation of this parameter allow to obtain the value corresponding to ntf = 3 or ntf = 4.

We prefer using ntf = 4 that gives lower noise, no overshoot on the time-surface temperature curve and reasonable heat flux density. To improve the results of the inverse method, some modifications of the instrumentation would be necessary, but present technical difficulties: the decrease of the sampling period, dt, is difficult in the context of unattended acquisition; it may require the diminution of the number of channels. The decrease of the location's depth of the first thermocouple is fastidious due to the proximity of the surface and of the tolerance of the drilling operation.

The time response of the thermocouples, τ ThC , is also a parameter that could generate some inaccuracy in the results of the inverse method. The response of the inverse method has been tested as a function of a dimensionless parameter τ ThC /τ . The sensibility study [START_REF] Dour | Development of a nonintrusive heat transfer coefficient gauge and its application to high pressure die casting: Effect of the process parameters[END_REF] has shown that the maximal value of the heat flux density is reached if this parameter is lower than 10%. Considering the thermocouple as a first-order system, its response to a heat flux step function has then an origin slope equal to the ratio q max /τ ThC (q max is supposed to be the set point). Assuming that the origin slope is the maximal slope of the time-heat flux density curve, qmax , it could be estimated a value of the time response of the thermocouple lower than The deterministic uncertainties on heat flux evaluation depends on temperature measurements, thermal properties of the mould and uncertainty in the location of the first thermocouple [START_REF] Dour | Recommendations and guidelines for the performance of accurate heat transfer measurements in rapid forming processes[END_REF]. A normalized analysis of the thermal problem enables the direct estimation of the deterministic errors, as discussed in reference [START_REF] Dour | Recommendations and guidelines for the performance of accurate heat transfer measurements in rapid forming processes[END_REF], by derivation of the equations of the direct normalized problem. The error in the heat flux density dq could be determined from Eq. [START_REF] Höhne | Measurement and mathematical modeling of the heat transfer in the glass forming process, in consideration of the heat transfer coefficients and radiation influences[END_REF].

dq q = dT T + dk k + 1 2 dκ κ + q max k T dz ( 5
)
where T is the temperature range (taken to 150 • C according to Fig. 8), dT is the error in the temperature measurement (±1.5 • C according to the IEC 584 norm for Class 1 K-type thermocouples), dk (±1 W m -1 K -1 ) is the error in the thermal conductivity value assumed to be constant at 26 W m -1 K -1 , dκ (±0.5 × 10 -6 m 2 s -1 ) is the error in the thermal diffusivity value assumed to be constant at 4.9 × 10 -6 m 2 s -1 , q max is the peak value of the evaluated heat flux density (3.2 MW m -2 ) and dz is the inaccuracy in the location of the thermocouple estimated to ±0.055 mm, according to the tolerance in the drilling operation (±0.03 mm) and the accuracy of the location of the thermocouple (0.25 mm diameter) in a greater hole (0.3 mm diameter). With these values, the deterministic inaccuracy in the heat flux density estimation was found to be ±14.5%. The most influent factors are the diffusivity (5.1%), the location of the thermocouple (4.5%) and the conductivity (3.8%) of the steel, and finally the temperature measurement (1%). The verifications made above concerning the ntf parameter and the time response of the thermocouples does not allow to validate the results of the inverse method. For that use, the measurement of the thermocouple located at 7 mm from the surface of the gauge is compared to the temperature calculated by using the evaluated heat flux density as boundary conditions. The comparison of the temperatures measured and calculated at 13 and 1 mm from the surface are not discussed because they are close to zero, the first one being imposed and the second one serving as optimization objective. At 7 mm, the error is lower than 20 • C at each gauge (see Fig. 8, Section 3.3), what corresponds to a relative error about 10%. This error is still acceptable because it maximizes the deviation between experiment and calculation over the whole temperature field.

The reasons of such a deviation are of two types. First of all, the assumption is made in the thermal quadrupoles method that the heat transfer is one-dimensional and perpendicular to the surface. This hypothesis does not seem to be of great influence because the same deviation is observable at each gauge, whereas it should be more valid at gauge C1 where heat transfer parallel to the surface should occur considering the symmetry, than at gauge C2 and C3. To compare the velocity of the front matter and the diffusion of the heat in the bulk of the mould, the Peclet number would be useful. From a previous work [START_REF] Dusserre | Thermo-mechanical stresses in cast steel dies during glass pressing process[END_REF], the velocity of the front matter when the glass flow on the plate surface of the mould has been estimated about 0.25 m s -1 . The thickness of the gob is then close to 0.03 m which is also close to the thickness of the mould. Using these values and a diffusivity of 4.9 × 10 -6 m 2 s -1 , one obtains a value of the Peclet number about 1530. With such a high value, the velocity of the front matter is widely higher as the diffusion of the heat and the heat flux density could be supposed to be applied instantaneously on the surface of the mould.

The second reason which could explain those deviations is the assumption that the material properties are independent of the temperature. Fig. 6 shows the differences between the thermal diffusivity and conductivity of the GX30Cr13 steel and their mean value considered in the range of use of the inverse method. The error made on the thermal conductivity is very weak (lower than 4%), but an error of about 10% exists on the thermal diffusivity for the highest and the lowest temperatures of the range considered.

Choice of a representative cycle

The repeatability of the thermal cycles studied by the described method is quantified using the maximal and minimal values evaluated for 31 stabilized cycles (see Fig. 7). The peak mean values are respectively 2.2 MW m -2 during pressing and -0.220 MW m -2 during mould cooling. The standard deviation of the maximal flux value is nevertheless only twice the one of the minimal value whereas the mean absolute values are about ten times higher (see Table 3). This is due to the higher sensitivity to the noise for the standard deviation of the weaker values. From this statistical study of the repeatability of the consecutive cycles, a cycle representative of a standard production process is chosen.

Upper side temperature surface and heat flux density

Fig. 8 shows the temperature measurements at the three gauges for the typical average cycle. On the right side of the figure set is shown the heat flux density evaluated for the three gauges. All the key results (T min , T max , T , q) are summarized in Table 4.

The maximal temperature at the surface of the mould varies between 565 and 605 • C, depending on the location of the gauge (see Table 4). It is much lower at gauge C1 (middle of the mould) than at the gauge C2 and C3, which could be explained by a more efficient cooling due to the centering of the upper and lower cooling setups (the minimal temperature is indeed lower as well). No significant differences are observable on the temperature amplitude between the minimum and the maximum, but it is slightly lower at gauge C3 (120 • C instead of about 140 • C at gauges C1 and C2). The small difference may be attributed to the lower heat absorbed at this location due to the shorter contact time (see Fig. 8). The temperature amplitude at the lower side of the mould is about 25 • C, between 335 and 360 • C.

At the delivery station (n • 1), the temperature increases very quickly at gauges C1 and C2, but not at gauge C3 because the latter is not located under the gob. The cooling of the glass generates a thermal contraction having for consequence that the gob leans on its own periphery. The contact is then effective at gauge C2 and non-effective at gauge C1, explaining the higher value of the heat flux density at gauge C2 (3.2 instead of 2.2 MW m -2 ).

While pressing the glass, the contact pressure increases and a second peak of flux appears at gauge C1 (0.9 MW m -2 ). This second peak does not exist at gauge C2 because the contact was still fully effective before pressing. A first peak is observable at C3 when the glass flows in front of this gauge. The value of this peak is about the same as at gauge C2 during the waiting time (3.0 MW m -2 ) but it is not as wide because of a shorter contact time.

At the end of the pressing stage, the temperature is maximal at each gauge. Once the upper mould is moved up, the contact pressure is much weaker and the surface temperature begins to decrease. At this time, the behavior of the gauge C3 is sensibly different than the other two. A very quick surface temperature diminution of about 30 • C takes place followed by a slower increase of about 20 • C, whereas the temperature decreases almost regularly at gauges C1 and C2. A possible reason of such a behavior is a detachment between the glass and the mould, more important at the periphery than at the center due to a deformation of the mould when the upper mould moves up. After that, the contact pressure and then the temperature could increase due to a creep of the glass under its own weight or due to the upper air blast. The modeling of this phenomenon would be very difficult because of the number of parameters controlling it, such as glass mechanical behavior, temperature field in the glass, behavior of the contact between glass and mould . . . The surface temperature of both upper and lower mould could play an important role on the deformation of the glass when the upper tool moves up. Indeed, the separation force between a glass melt and a metallic substrate increases exponentially for mould temperatures above 550 • C whatever the mould material, as shown by Manns et al. [START_REF] Manns | Glass in contact with mould materials for container production[END_REF].

Comparison between the upper and lower cooling

Fig. 8 shows that the negative heat flux density oscillates nearby two values. The first value occurring during the stages n • 1 (before delivery), 10 (after vacuum take out), 11, 12 and 16 corresponds to a free convection cooling. The second value, lower than the first one, occurs during the stages n • 13, 14 and 15 and corresponds to forced convection while upper air blast cooling.

The mean cooling heat flux density on the upper side of the mould is about -50 kW m -2 for the free convection (a very weak difference exists between each gauge, see Table 4) and between -137 kW m -2 and -112 kW m -2 for the forced convection (the value is higher at the center of the mould due to the centering of the cooling setup). An important standard deviation of about 23 kW m -2 is observed on the mean values above, due to the non-negligible noise existing in the results of the inverse method. From these heat flux density one can estimate the convection coefficient between the mould temperature

Table 4

Characteristics of a representative thermal cycle at gauges C1, C2 and C3 and the air at 35 • C. One would find respectively 118, 97 and 89 W m -2 K -1 for free convection at gauges C1, C2 and C3 (the mould temperature is respectively of 500, 520 and 540 • C) and 330, 285 and 241 W m -2 K -1 for forced convection (the mould temperature is respectively of 450, 480 and 500 • C).

T max T min T Peak of flux, MW m -2 q free q forced h up free h up forced • C • C • C S t a t i o nn • 1 S t a t i o nn • 2 k W m -2 kW m -2 W m -2 K -1 W m -2 K -1 C1
Q t 1 →t 2 = t 2 t 1 q(t) dt (6)
The heat density Q + (defined by Eq. ( 6)) delivered by the gob to the mould during one stabilized cycle (see Table 5) is larger than the heat density Q -absorbed by convection on the upper side of the mould during the same stabilized cycle. This demonstrates the importance of the lower side cooling, twice higher as the upper side one. This could be explained by the duration of the cooling (during the whole cycle on the lower side, but only between vacuum take out and delivery on the upper side), and the presence of cooling fins on the lower side of the mould. The convection coefficient on the lower side could be estimated considering the heat quantity Q + + Q -to be absorbed by convection at the lower side during a cycle (48 s) to obtain an equilibrated heat balance. Considering a lower side mould temperature of about 310 • C and an ambient temperature of 35 • C, the values estimated are respectively of 332, 408 and 320 W m -2 K -1 at gauges C1, C2 and C3. These values are sensibly higher to the ones estimated for forced convection on the upper side, but this difference is understandable due to the presence of cooling fins.

Application to finite element calculation

Temperature field in the whole mould

The experimental results described above are of great interest for the calculation of the thermo-mechanical loading in- side the whole bulk of the mould. A weak coupled thermomechanical calculation has been realized with the finite element software ABAQUS. In a first approach, the mould was supposed to be axis-symmetric with a projected working surface equivalent to the elliptical geometry. The meshing was constituted of 4834 quadrangular elements and 5404 nodes (linear interpolation of the temperature), see Fig. 9. The uniform initial condition corresponds to the preheating temperature (370 • C).

Q + Q - Q + + Q - h low forced h abaqus MJ m -2 MJ m -2 MJ m -2 W m -2 K -1 W m -2 K -1 C1
The boundary conditions applied on the meshing are those estimated by inverse method from the temperature measurement in production conditions. Two kinds of boundary conditions are applied: heat flux on the working surface and convection coefficient on the external sides.

The work surface has been divided into 17 faces. On the first one, the heat flux density estimated at gauge C1 is applied. On the second one, the heat flux density estimated at gauge C2 is applied. On the third face, the heat flux density estimated at gauge C2 is applied with a time lag in order to simulate the creeping of the gob under its own weight before pressing. On the fourth face, the heat flux density estimated at gauge C3 is applied. On the faces 5 to 17, the heat flux density estimated at gauge C3 is applied with a time lag in order to take into account the progression of the glass front while forming. The size of the different faces and the value of each time lag have been determined from the results of a simulation with FORGE2 FEM package of the glass flow during the glass pressing process. A time lag corresponds to a time step of this calculation, and the size of the corresponding face is equal to the progression of the front matter during this step. The initial and boundary conditions used in this previous glass flow simulation were taken from the literature (time dependent heat exchange coefficient according to the results of the reference [START_REF] Moreau | Méthodologie expérimentale pour la détermination du coefficient de transfert thermique à l'interface verre/outil[END_REF]) and does not match exactly our processing conditions. Nevertheless, the cinematic of the press was imposed, and even if the stress results may be discussable, the movement of the front matter should be well described. The details of this simulation are available in a previous paper [START_REF] Dusserre | Thermo-mechanical stresses in cast steel dies during glass pressing process[END_REF].

The external surface of the mould has been divided into 4 faces. The first three are located under each gauge and the convection coefficients applied have been estimated previously from temperature measurement (see Table 5). On the last face corresponding to the lateral external side of the mould, a convection coefficient of 110 W m -2 K -1 has been applied, as estimated above for free convection (see Table 4).

Sixty cycles have been simulated with these boundary conditions in order to obtain a stabilization of the mean temperature. A significant difference between the measured and calculated stabilized temperature has been observed due to the one dimensional assumption used to determine the values from the experimental data, the ratio between heat transfers parallel and perpendicular to the surface increasing with the distance of the surface. The values of the convection coefficient applied under each gauge have then been optimized in order to minimize this difference. The optimized values are respectively of 560, 180 and 160 W m -2 K -1 instead of respectively 332, 408 and 320 W m -2 K -1 . After optimization, the differences between experience and simulation are very weak at each comparison points (measured bulk temperatures, estimated upper surface temperature and measured lower surface temperature). As shown in Fig. 10, the differences between measurements and simulation are very weak (less than 20 • C at each comparison point over the whole cycle). The surface temperature is very well evaluated. The temperature dependence of the thermal properties of the steel has been taken into account in the finite element calculation leading to a better agreement with the experimental data at 7 mm with finite element (Fig. 10), than with thermal quadrupole (Fig. 8).

Validation of the simulation of the glass forming by pressing

A simulation of the pressing process (see [START_REF] Dusserre | Thermo-mechanical stresses in cast steel dies during glass pressing process[END_REF] for details) has been realized with the finite element software FORGE2. The simplified thermal boundary conditions applied on this fully coupled thermo-mechanical calculation in the gob and in the lower mould were a heat exchange coefficient and eventually a time dependent opposite surface temperature (the boundary conditions described further could differ from those presented in [START_REF] Dusserre | Thermo-mechanical stresses in cast steel dies during glass pressing process[END_REF], due to the availability of experimental data).

During waiting time before pressing, the heat exchange coefficient and the opposite temperature (corresponding to gob surface temperature) applied are function of time as described in Fig. 11. The temperature of the gob surface is computed independently by imposing the heat flux density estimated by inverse method on the glass with a one-dimensional finite differences calculation (the heat exchange coefficient can then be estimated). During the pressing stage and the glass cooling stages, the heat transfer coefficient is supposed to be constant (7 kWm -1 K -1 for pressing and 500 W m -1 K -1 for glass cooling) and the opposite temperature is the one calculated at the surface of the gob mesh. The initial condition results from the independent gob temperature simulation during the waiting time. During mould cooling, the values of the heat exchange coefficient and opposite temperature are those estimated at gauge C1 (118 W m -1 K -1 for free convection, 330 W m -1 K -1 for forced convection and ambient temperature of 35 • C).

The results presented in Fig. 11 compare the surface temperature calculated with FORGE2 and the surface temperature estimated at gauge C1 with the inverse method. It can be observed a good description of the evolution of the temperature during waiting time and pressing, as well as while cooling the mould surface by free or forced convection, even if a significant over evaluation exists. During the cooling of the glass part, the evolution of the temperature is not realistically described by the hypothesis of constant heat transfer coefficient explaining partially the later over evaluation of the surface temperature. This error may be due to the very complex behavior of the contact between the glass and the mould after pressing as discussed at the end of Section 3.3.

Conclusion

The results of an investigation of the heat flux density extracted from a gob by glass pressing in a mould made out of martensitic stainless steel was presented. With this aim in view, an instrumentation as non-intrusive as possible has been developed to measure temperatures in the bulk of the mould in normal production conditions. This instrumentation is composed of three gauges, each allowing temperature measurements at 1, 7 and 13 mm from the surface of the mould. An inverse method based on the thermal quadrupoles and the Beck's optimization method has been used to exploit these measurements and evaluate the heat flux density absorbed as well as the temperature at the surface of the mould. The study of the difference between each gauge during a representative cycle shows the influence of the production stages, the heterogeneity and the complexity of the thermal loadings at the surface of the mould. The pre-dominance of the lower air blast cooling over the upper one has been demonstrated comparing the heat quantities delivered by the gob and absorbed by convection on the upper side of the mould. This observation justifies the presence of cooling fins on the lower surface. The quality of these results has been discussed in terms of the validity of the hypothesis and of the influence of the intrinsic parameters of the inverse method and of the instrumentation. The use of the heat flux density evaluated as boundary conditions in a finite element calculation allow to calculate with a very good accuracy the evolution of the temperature field in the whole bulk of the mould. It has also been shown that the hypothesis of constant heat transfer coefficient is not sufficient to describe properly the evolution of the surface temperature of the mould, especially during the glass part cooling stages. This also reflects the complexity of the glass/mould interface as deduced from our measurements.

CFig. 1 .

 1 Fig. 1. Different stages of the glass pressing process: upper and lower cooling air flow.

Fig. 2 .

 2 Fig. 2. The instrumented mould in its industrial context just before the delivery.

Fig. 3 .

 3 Fig. 3. Description of the instrumented mould: location of the gauges and thermocouples housings.

Fig. 4 .

 4 Fig. 4. Relevance domain of the inverse method [25,27].

Fig. 5 .

 5 Fig. 5. Influence of the ntf parameters from 2 to 5 (max heat flux density and noise).

  cycle and the duration of the heat flux application, supposed to follow a step function. Using the characteristics of the cycle used to show the influence of the ntf parameter (variation of the surface temperature at gauge C3 about 110 • C, duration of the heat flux density step function of 0.12 s, thickness of the mould and properties of GX30Cr13 steel), one obtains a value about 3.9 MW m -2

Fig. 6 .

 6 Fig. 6. Evolution of the thermal diffusivity and conductivity of the GX30Cr13 steel with temperature and mean value used for inverse method (Aubert & Duval data).

Fig. 7 .

 7 Fig. 7. Successive heat flux density peaks in the production steady state at gauge C1.

Fig. 8 .

 8 Fig. 8. Results of the temperature measurements (on left hand) and the inverse method determination of the heat flux density (on right hand) in Sunderland at gauges C1 (top), C2 (middle) and C3 (down). The temperatures computed with thermal quadrupoles (line) and measured (dot) at 7 mm are compared. The vertical rules materialize the transition between two consecutive stations, as numbered on Fig. 1.

Fig. 9 .

 9 Fig. 9. Meshing of the mould for finite element calculation (Abaqus).

Fig. 10 .

 10 Fig. 10. Comparison between the measured temperature (on left hand) and the temperatures computed with finite element method (on right hand) in Sunderland at gauges C1 (top), C2 (middle) and C3 (down). The vertical rules materialize the transition between two consecutive stations, as numbered on Fig. 1.

Fig. 11 .

 11 Fig. 11. FORGE2 calculation: boundary conditions during waiting time before pressing (on left hand) and comparison of the surface temperature calculated and estimated by inverse method (on right hand).

Table 1

 1 Chemical composition (wt.%) of GX30Cr13 steel (designation NF EN 10027-1)

	C	C r	S i	M n	N i	F e
	0.3	13	0.6	0.9	<0.5	balance

Table 2

 2 Elaboration, heat treatment and hardness of GX30Cr13 steel

	Sand casting	Normalization	Austenitizing	Quenching	Tempering	Hardness
	1600 • C/Ar	880 • C/3 h	1000 • C/8 h	Air blast	680 • C/6 h	280 HB/27 HRc

Table 3

 3 Repeatability of the thermal cycles at gauge C1 (population: 31 cycles)

		Mean	Standard	Maximal	Minimal
		value	deviation	value	value
	q max , MW m -2 q min , MW m -2	2.20 -0.220	0.0442 0.0213	2.31 -0.178	2.11 -0.250

Table 5

 5 Lower side convection coefficient estimated at gauges C1, C2 and C3 (T mould = 310 • C; T air = 35 • C; t cycle = 48 s)

Acknowledgements

The authors greatly acknowledge the financial support of this work by Arc International Cookware, Aubert & Duval and Saint-Gobain SEVA. This study would not have been carried out without the participation of the production teams of the manufacturing plants of Châteauroux and Sunderland. Special thanks are expressed to M. Daniel Deletang.