N Gascoin

L Romagnosi

G Fau

P Gillard

Pyrolysis in Porous Media: Numerical Analysis and Comparison to Experiments

Only limited studies are available experimentally to investigate hydrocarbon fuel pyrolysis, which can be of practical interest for the active cooling of head loaded components in aerospace vehicles such as combustors in rocket engines. The numerical simulation is an additional way to investigate the related phenomena (heat and mass transfer with chemistry). After a first code validation with experiments based upon inert gases, the code is further extended towards permeation with reactive dodecane and a reasonable agreement is found. The experimental accuracy of the permeability's' determination is confirmed numerically to be within ± 30 %. Numerically it is shown that this accuracy is due to strong flow spatial heterogeneities. The border effect of the test cell is found to be related to the permeable medium thickness whereas the temperature field is correlated to the reaction zone. Two different kinetic mechanisms are used to investigate the chemistry effect on the heat and mass transfer. They provide also a better analysis of the fuel pyrolysis and of the products' formation.

I. Introduction

Studying porous flow is of interest in numerous fields of research and more specifically in the present framework of active cooling intended to be used on highly thermally loaded structures [START_REF] Langener | Experimental Investigation of Transpiration Cooling with Subsonic and Supersonic Flows at Moderate Temperature Levels[END_REF][START_REF] Kuhn | Application of Transpiration Cooling for Hot Structures, RESPACE -Key Tech. for Reusable Space Sys[END_REF]. Transpiration cooling through a porous structure is highly efficient as shown by Langener et al. [START_REF] Langener | Experimental Investigations on transpiration Cooling for Scramjet Applications Using Different Coolants[END_REF][START_REF] Cheuret | Numerical Investigation on Transpiration Cooling for Scramjet Applications using Different Coolants[END_REF]. This technology is found in rocket engines and chambers for example [5,6]. This reduces the temperature of the solid permeable materials used to sustain the large heat load encountered by the combustor [START_REF] Fau | Fuel Pyrolysis through Porous Media: Coke Formation and Coupled effect on Permeability[END_REF][START_REF] Gascoin | Supersonic combustion of hydrocarbons pyrolysed mixture with detailed chemistry[END_REF]. Carbon-or oxide-based composite materials are preferred for thermal, mechanical and chemical reasons [START_REF] Bouchez | PTAH-SOCAR Fuel-Cooled Composite Materials Structure: 2009 Status, 16th AIAA/DLR/DGLR International Space Planes and Hypersonic Systems and Technologies Conference[END_REF]. In addition, when using the fuel as a coolant, the so-called regenerative technique, a pyrolysis is generally encountered for fluid temperature over 800 K [START_REF] Gascoin | Etude et mesure de paramètres pertinents dans un écoulement réactif[END_REF]. This generates endothermic reactions which contribute to the heat flux absorption. The formation of a thin film inside the combustion chamber, the so-called film cooling technique, is another parameter participating in the structure cooling [START_REF] Arnold | Film Cooling of Accelerated Flow in a Subscale Combustion Chamber[END_REF]. The main drawbacks of these coupled solutions are the multiphysics' interactions which modify transiently the permeabilities of the structure and the cooling efficiency [START_REF] Fau | Fuel Pyrolysis through Porous Media: Coke Formation and Coupled effect on Permeability[END_REF]. The fuel cracking produces light and heavy species which are coke precursors [START_REF] Gascoin | Characterisation of coking activity during supercritical hydrocarbon pyrolysis[END_REF]. The formation of solid particles modifies the fluid filtration and consequently the heat and mass transfers [START_REF] Gascoin | Thermal and hydraulic effects of coke deposit in hydrocarbon pyrolysis process[END_REF].

Experimental chemical studies generally focus on the fuel pyrolysis as a function of the fuel nature, of the operating conditions or of the reactor's nature [START_REF] Pant | Pyrolysis of n-heptane : kinetics and modeling[END_REF][START_REF] Zeppieri | Modeling concepts for larger carbon number alkanes : a partially reduced skeletal mechanism for n-decane oxidation and pyrolysis[END_REF][START_REF] Herbinet | Thermal decomposition of n-dodecane : Experiments and kinetic modeling[END_REF][START_REF] Yu | Supercritical-phase thermal decomposition of binary mixtures of jet fuel model compounds[END_REF][START_REF] Edwards | Fuel Composition Influence on Deposition in Endothermic Fuels[END_REF]. Some numerical works are also available [START_REF] Gascoin | Real-time method for the identification and quantification of hydrocarbon pyrolysis products: Part II. Application to transient pyrolysis and validation by numerical simulation[END_REF][START_REF] Slavinskaya | PAH and Soot Formation within Kerosene Pyrolysis[END_REF]. Other studies are related to the heat and mass transfers involved in the active cooling by means of computing [START_REF] Lina | Experimental and Numerical Investigation of Active Cooling Ceramic Matrix Composite for Ramjet Propulsion System[END_REF] and of testing [START_REF] Langener | Experimental Investigation of Transpiration Cooling with Subsonic and Supersonic Flows at Moderate Temperature Levels[END_REF][START_REF] Kuhn | Application of Transpiration Cooling for Hot Structures, RESPACE -Key Tech. for Reusable Space Sys[END_REF]. Nevertheless, there is only rare studies aiming at grouping both experimental and numerical approach to provide a deep analysis of the multiphysics' phenomena (thermal, fluid, mechanical and chemical points) [START_REF] Fau | Fuel Pyrolysis through Porous Media: Coke Formation and Coupled effect on Permeability[END_REF]. The complementarity of these two kinds of work is due to the fact that numerical studies still face strong difficulties in solving massively detailed kinetic mechanism and in handling porous flow with microscopic description over full length scales and over long timescales. On the opposite, it is difficult experimentally to determine microscopic or local parameters such as the temperature or the species concentration and to measure the heat flux or the residence time. For this reason, a European collaboration has been settled to establish a numerical modelling of reactive porous flow with semi-detailed chemistry and taking the other physical parameters into account (heat and mass transfers, multi-diffusion, porous flow modelling, compressibility effects, real fluid properties). The authors here aim at comparing the numerical results to the experimental ones on dodecane pyrolysis cases in order to determine the model sensitivity and accuracy. Numerical data are of major importance for the analysis of the experimental results. Dodecane is used to compare with experiments but methane should be used in latter steps for rocket combustors application.

II. Material and Methods

A. Numerical simulation

A 2-D geometry with 14,242 structured cells (Figure 1) models the test cell to be detailed in next section. The solver handles compressible steady-state configurations. The Navier-Stokes equations (Eq. 1 to Eq. 3) are modified taking into account the porosity effect of the sample, the consequent pressure drop through the medium and the effective thermal conductivity of the porous zone (fluid plus porous material). The Forchheimer's equation is used to simulate the pressure drop via a sink term in the momentum equation (Eq.4), where α and C 2 are respectively 2.061•10 -13 m 2 and 4924000 m -1 (obtained experimentally [START_REF] You | Detailed and simplified kinetic models of n-dodecane oxidation: The role of fuel cracking in aliphatic hydrocarbon combustion[END_REF]).

   

  i f f f S g P v v t v                        (1)     0       v t f f    (2)               h f i i i eff f f s s f f S v J h T k P E v E E t                          1 (3) i i i v v C v S    2 1 2    (4)
The pressure-based solver guaranties a faster and stable convergence. The boundary conditions correspond to constant: inlet mass flow rate, pressure outlet and wall temperature. The conditions measured experimentally are used directly in the software to simulate identical test cases. The simulation settings change depending on the case under study: nature of the fluid (nitrogen or kerosene surrogate), operating temperatures (from 300 K to 1200 K) and mass flow rate (60 mg.s -1 to 3.98 g.s -1 ). Additional details will be provided when necessary for each test case (section III).

Figure 1 Computational domain layout.

A detailed kinetic mechanism (56 species and 289 reactions) simulates the dodecane pyrolysis [START_REF] You | Detailed and simplified kinetic models of n-dodecane oxidation: The role of fuel cracking in aliphatic hydrocarbon combustion[END_REF]. The present authors have reduced this mechanism (15 species, 20 reactions) for computation cost reason with limited loss of accuracy. Two primary reactions ("1" and "2" in Table 1) allow the cracking of dodecane producing ethylene (C 2 H 4 ) and two radicals but-1-yl (pC 4 H 9 ) and n-propyl (nC 3 H 7 ). These two reactions controlling the pyrolysis process correspond to a high activation energy (≈ 3.68•10 5 J.mol -1 ). The sensitivity analysis [START_REF] You | Detailed and simplified kinetic models of n-dodecane oxidation: The role of fuel cracking in aliphatic hydrocarbon combustion[END_REF] for n-dodecane pyrolysis shows that the fuel conversion rate is controlled by C-C fission of n-dodecane and H-abstractions by H and CH 3 . So the present authors primary considered adding at least the reactions allowing the formation of these two radicals (monatomic hydrogen via backward reaction "4" and "5" and methyl radical via "12"). The presence of H and CH 3 accelerates the pyrolysis rate via propagation reactions for reactions "6" → "10" and "13" → "17". The ramification is indirect via dissociation of light species as hydrogen H 2 (reaction "11") and ethane C 2 H 6 (reaction "20"). The reactions involving directly the interaction between the n-dodecane and the radicals H and CH 3 include both Habstraction and β-scission of n-C12-alkyl-radical (C 12 H 25 ).

By comparing this reduced mechanism under perfectly stirred reactor configurations at 3.5 MPa with the original one [START_REF] You | Detailed and simplified kinetic models of n-dodecane oxidation: The role of fuel cracking in aliphatic hydrocarbon combustion[END_REF] and with a very detailed one (271 species and 1449 reactions) [23], the present reduced scheme is judged adequate for the specific conditions used in this paper despite important discrepancies (Figure 2). The overall performance is generally poor but they are acceptable in the specific thermal range 900 K-1200 K at 35 bar and for residence time comparable to those in the system. At 900 K and after 0.7 s of simulation (which is roughly the residence time found later in the cell, see section C), the conversion of dodecane is adequately simulated by the reduced scheme (Figure 2a) while the initial detailed one of You et al. [START_REF] You | Detailed and simplified kinetic models of n-dodecane oxidation: The role of fuel cracking in aliphatic hydrocarbon combustion[END_REF] presents a factor 2 compared to the accurate reference mechanism of Herbinet et al. [23]. A factor 2 is found on ethylene between the reference and the reduced scheme, which shows that the products distribution is predicting with a limited accuracy. After 10 s at 900 K (constant pressure of 3.5 MPa), ethylene and ethane are reproduced fairly contrary to some species like hexene (Figure 2b). At 1200 K and 35 bar, after 10 s, the discrepancies are still important between the reference and the reduced mechanisms, mainly concerning the methane formation while the original mechanism of You et al. [START_REF] You | Detailed and simplified kinetic models of n-dodecane oxidation: The role of fuel cracking in aliphatic hydrocarbon combustion[END_REF] gives satisfactory results. The major compound found in the reduced mechanism is ethylene (Figure 2c). Looking at the dynamics of the pyrolysis (Figure 2d), the original mechanism presents a shift on the log-time scale compared to the reference mechanism. The reduced mechanism presents a mean behaviour around the reference one. As a consequence, in order to perform first kinetic computations within 2-D geometry, the reduced mechanism is used but one should pay attention to the limited accuracy in terms of products formation. This mechanism has been optimised for representing the dodecane consumption and it cannot accurately represent the formation of radicals, of light species like hydrogen and of middle-weight species like hexene. This reduced mechanism is a compromise which should be later improved by increasing the number of species and of reactions and, as a consequence, the computation cost.

Table 1 Reduced kinetic mechanism used in this paper. [1]

NC12H26 => 3C2H4+2nC3H7

[12] 2CH3<=>H+C2H5 [2] NC12H26 => 2C2H4+2pC4H9 [13] NC12H26+CH3=>4C2H4+pC4H9+CH4 [3] C2H4+C2H5<=>pC4H9 [14] NC12H26+CH3=>C4H81+2C2H4+pC4H9+CH4 [4] C6H12+H = C2H4+pC4H9 [15] NC12H26+CH3=>C3H6+C6H12+nC3H7+CH4 [5] C5H10+H = C2H4+nC3H7 [16] NC12H26+CH3=>C5H10+2C2H4+nC3H7+CH4 [6] NC12H26+H => 4C2H4+pC4H9+H2 [17] NC12H26+CH3=>C6H12+C2H4+pC4H9+CH4 [7] NC12H26+H => C4H81+2C2H4+pC4H9+H2 [18] C2H6+CH3<=>C2H5+CH4 [8] NC12H26+H => C3H6+C6H12+nC3H7+H2 [19] CH3+CH2<=>C2H4+H [9] NC12H26+H => C5H10+2C2H4+nC3H7+H2 [20] 2CH3(+M)<=>C2H6(+M) [10] NC12H26+H => C6H12+C2H4+pC4H9+H2 LOW / 1.

B. Experimental setup

The COMPARER pyrolysis test bench is used to pressurize and to heat the fuel under flow conditions [START_REF] Gascoin | Test Bench Dimensioned by Specific Numerical Tool[END_REF]. Its main characteristics are the following:

 Maximum operating conditions: 1800 K, 80 bar, 0.6 g.s -1 for liquid fluid and 6 g.s -1 for gas  Chemical analysis: Gas Chromatograph/Mass Spectrometer-Flame Ionisation Detector-Thermal Conductivity Detector with10 sampling loops for a posteriori quantitative analysis, Fourier Transform Infra Red spectrometer (inline direct quantitative analysis)  Sensors: 5 mass flow rates, 5 pressure transducers, over 10 K-type and R-type thermocouples with data acquisition system (16 bits, 48 channels, 0.1 Hz) A permeation test cell contains the porous sample (Figure 3). This cell is inserted inside the furnace of the COMPARER bench and it is connected to the fluid supply system and to the appropriate sensors. The permeable medium bounds the cell in two high and low pressure chambers (upstream and downstream to the porous medium respectively). These chambers are noted HPC and LPC in this paper. An inlet pipe provides the fuel into the system. This cell is connected to a dynamic sampling system to get hot pressurized samples at three location points in the cell. These samples are later cooled and chemically quantified by the GC/MS-FID-TCD. This allows having the longitudinal pyrolysis profile of the fuel under given operating conditions. The FTIR analysis provides real-time quantification at the process outlet on the gas products only after having cooled and expanded the pyrolysis mixture at ambient conditions.

The permeable medium can be changed. In the present work, an isotropic stainless steel material is preferred to composite one to avoid considering complex microstructure (fibres, layers). It is characterized by a porosity around 30 %, a grain diameter of 14.1 µm and a pore diameter of 4.1 µm [START_REF] Gascoin | High temperature and pressure reactive flows through porous media[END_REF]. Further geometrical information can be found in [START_REF] Gascoin | High temperature and pressure reactive flows through porous media[END_REF].

III. Fuel pyrolysis analysis A. Preliminary validation on inert test cases and first experiments analysis

Nitrogen is injected with constant mass flow rate and the inlet pressure is 2.5 MPa. The validation is made on the pressure drop through the porous material. Several temperatures have been tested from 300 K to 600 K. The thermal conductivity of the stainless steel -nitrogen system (k eff ) is evaluated as a function of temperature using an averaged method between Parallel model and Maxwell Upper Limit [START_REF] Travajno | Calcul multiphysique dans une paroi poreuse[END_REF].

The simulations performed with the initial C 2 value presented some discrepancies compared to the experiments. Thus, the authors have corrected the C 2 value by a multiplication factor of 1.25 (new C 2 = 6155000 m -1 ) because the Forchheimer accuracy is ± 30 % [START_REF] Gascoin | High temperature and pressure reactive flows through porous media[END_REF]. As a result, the computations are in good agreement with the experimental results (Figure 4). Note that, since the Forchheimer's coefficient is linked to the square of the velocity, the changes in C 2 have a negligible effect on the pressure drop for low mass flow rates.

The numerical simulation interest is clear because it can help validating the permeabilities values and understanding the reasons of any possible error. For example here, the complex velocity and pressure field upstream to the porous medium (Figure 5) generates a recirculation zone, which makes the pressure measurements sensitive to the spatial position. This impacts the experimental estimation of K D and K F . As a consequence, it is wanted to verify how the pressure measures are appropriate for the kerosene experiments. A similar cold simulation has thus been done using kerosene surrogate (n-dodecane), which is liquid at ambient temperature (T inlet = 300 K). The mass flow rate is 65 mg.s -1 . The density is fixed at 746 kg.m -3 and the inlet pressure is 3.5 MPa. The high viscosity of the kerosene (µ = 1.37•10 -3 Pa.s) and its low injection velocity bounds the recirculation zones at the cell inlet between x ≈ -0.025 m and x ≈ -0.018 m (Figure 6). This means that the flow invests the porous medium almost uniformly and, in this case, it is not so important where the pressure sensors are located.

B. Validation with reactive dodecane experiment

The experiments were achieved by Fau et al. [START_REF] Fau | Fuel Pyrolysis through Porous Media: Coke Formation and Coupled effect on Permeability[END_REF] at 3.5 MPa with a varying mass flow rate between 60 mg.s -1 and 65 mg.s -1 due to stability problems. The setup furnace temperature was 900 K and the temperature experimentally measured at the outside wall of the cell was 890 K. The fluid temperature measured inside the test cell (in the high pressure chamber) was 820 K. Longitudinal thermal gradients measured experimentally were lower than 10 K roughly between two positions: 3 mm before the porous medium and 3 mm after, following the main fluid flow direction. The experiments (Figure 7) show a high dodecane content upstream to the porous medium (in the High Pressure Chamber) and inside the porous medium (83.3 wt.% and 81 wt.% respectively) while a mass fraction around 52 wt.% has been found downstream, in the Low Pressure Chamber. The experimental quantification of the main species considered by the reduced mechanism (Table 1) are summarized in Figure 7 for comparison purpose.

Numerically, it is difficult to represent all the experimental fluctuations. As a consequence, several test cases have been computed (Table 2). For all of them, the pressure is 3.5 MPa, the fluid is the dodecane and the reduced kinetic mechanism presented above is used (15 species, 20 reactions). The real gas model enables computing the fluid properties, which requires knowing the species critical properties such as the critical temperature, pressure and volume and the acentric factor. For the first case (Table 2), minimum (Figure 8a) and maximum (Figure 8b) mass fractions at several longitudinal x-axis positions are given for a selection of species. The upstream zone presents a low dodecane content (Figure 8a) due to the recirculation in which quasi static fluid is found. Hexene, pentene and ethylene are found in high proportion (Figure 8b). Comparing these results with the second case (Figure 8c,d) clearly shows that the higher temperature increases the dodecane pyrolysis because the minimum species quantities are much higher for this second case (Figure 8c) compared to the first one (Figure 8a). As a consequence, the maximum dodecane quantities are lower for this second case (Figure 8d). The same major products are found. The discrepancies between the two numerical cases are attributed in priority to the heat transfers but the mass flow rate also plays a role which will be investigated in section C. The dodecane mass fraction at the cell outlet of case 2 (x = 13 mm just before the outlet pipe) is found between 10 wt.% and 25 wt.%, which is lower than the experimental value. For the case 1, the extreme values are 15 wt.% and 85 wt.%. This includes the experimental value of 52 wt.%. Consequently, the test conditions are very important and slight changes (820 K to 890 K, in conjunction with mass flow rate variations) generate different results which are difficult to validate with the experiments. In addition, the kinetic mechanism has a role in the disagreement found between the numerical and the experimental data. The scheme has been tested (section II.A) and its accuracy is average. The pyrolysis process and the fluid dynamics are closely related because the density changes with the chemical composition and this affects the dynamics into the cell (residence time' and heat transfers' consequence). Hence, the kinetic mechanism becomes a critical parameter and adopting a detailed scheme is preferable.

Finally, because the case 1 is judged to be closer to the experimental conditions, the corresponding results are compared to the experimental ones (Figure 9). Three zones are defined: before the porous medium (for a position x = -3 mm), inside the porous zone (x = 1.5 mm) and at the process outlet (x = 30 mm). The mean rate is over-estimated in the HPC (Figure 9a) and in the porous zone (Figure 9b) by the computations because the spatial average over the y-axis for each x-axis position takes the stagnation zones near the wall into account while the experimental sampling locations are closer to the central axis of the cell. The agreement is much better at the bench outlet (Figure 9c). Nevertheless, due to the uncertainties presented above, there are too much parameters acting on the system to be certain that the numerical simulation is validated. This is particularly due to the kinetic scheme. The hexene content is overpredicted by the computations because the secondary mechanism of this species is not considered and it cannot be consumed. This explains the low quantities found for the light species like CH Experiments [START_REF] Langener | Experimental Investigations on transpiration Cooling for Scramjet Applications Using Different Coolants[END_REF] This numerical work c) Figure 9. Comparison of the pyrolysis products composition computed for case 1 and measured experimentally in the HPC (a), inside the porous zone (b) and at the process outlet (c). To verify the effect of the chemistry, a two-steps kinetic scheme is used to redo the same test case 2. A reduction of only 17 K of the fluid temperature between the inlet and the outlet is observed and the net heat flux at the wall is 3.2 W. This explains the very low pyrolysis rate at the outlet (remaining dodecane mass fraction of 98 wt. %). Consequently, the experimental values (pyrolysis rate of 48 wt. %) are bounded by the numerical ones obtained with two different kinetic schemes (pyrolysis rate of 80 wt. % for the detailed scheme and of 2 wt. % for the 2-steps scheme). The lower reactivity of the simplified mechanism is due to the radicals, such as methyl and monatomic hydrogen, which are not in the two-steps mechanism while they increase the reaction rate.

C. Use of the numerical simulation to investigate the reactive porous flow

Considering the third case (Table 2), the conversion rate of dodecane reaches 94 wt. % at the cell outlet (Figure 10a). A maximum thermal decrease of about 300 K is found (Figure 10b). This temperature fall is directly related to the enthalpy of reaction (Figure 10c), thus to the reaction zone. This point is particularly important because following the temperature would give precious indication on the chemistry. This is confirmed by looking at the molecular weight profile (Figure 10d). Such a result demonstrates the benefit of the computations compared to the experiments. This also shows that the recirculation zones, where the temperature is maximum (at y = 0.008 m and for x = -0.024 m and x = 0.013 mm), do not really participate to the pyrolysis. This is satisfactory when conducting pyrolysis in porous flow. This information concerning in which part of the system does the pyrolysis take place can be completed. Looking at the fuel consumption along the cell following the streamline on the axis of symmetry, a fast dodecane cracking occurs immediately at the inlet of the cell due the high inlet temperature (Figure 11a). The dodecane mass fraction entering the HPC is about 72 wt. % and it is around 66 wt. % at the porous medium inlet. It is 55 wt. % at its outlet, 28 wt. % at the LPC outlet and 18 wt. % at the pipe outlet. This means that the porous medium only represent 11 wt. % of the converted dodecane compared to the 82 wt. % which are converted in the entire system. Thus, the porous media represents about 14 % of the pyrolysis activity and its thickness (3 mm) corresponds to 4 % of the total length of the system (70 mm). In terms of residence time, the porous medium represents less than 10 %. The slope of the curve giving the residence time shows that the fuel is slowed down in the LPC because the slope increases. The fuel velocity increases in the outlet pipe due to the cross-section reduction.

Since the experimental mass flow rate slightly fluctuates between 60 mg.s -1 and 65 mg.s -1 , its effect has been investigated (case 3 and 4). The total residence time and its distribution are not so different (Figure 11b). Increasing the fuel mass flow rate decreases the pyrolysis rate at the outlet which means that the density and the residence time effect (acting similarly) have a lower impact than the convective heat transfers which increase (80 wt.% and 72 wt.% of pyrolysis rate for cases 3 and 4 respectively). This discrepancy of about 10 % is directly linked to the mass flow rate fluctuations of about 10 %. This demonstrates the direct effect of the fluid flow on the chemical parameters. Thus, simulating the experiments in which the flow rate varies is difficult numerically where the flow rate is fixed. This contributes to explain the disagreement found in section B. The effect of the temperature (comparison of cases 1 and 4) on the chemistry is represented by Figure 12. The strong chemical differences regarding the 2-D dodecane content (Figure 12a,b) are directly related to the thermal conditions. The residence time distribution is similar for case 1 (Figure 12c) than for cases 3 and 4 (Figure 11). Nevertheless, the residence time is about three times higher in case 1 than in other cases. This demonstrates that the temperature impacts the chemistry and the density, so the residence time, but its importance is much higher. A temperature decrease increases the density and decreases the velocity through the mass conservation. Thus, the residence time is higher but this is not balanced by an increased convective heat flux. This spate of events is important to analyse the results and this cannot be done experimentally, which shows the good complementarity of both tools.

D. Specific analysis of coupled phenomena in the porous medium

Comparing the dodecane conversion at the inlet and at the outlet of the porous medium for the fourth case, a coupling between the fluid flow and the chemistry is found (Figure 13a). Close to the symmetric axis (y = 0 mm), the flow speed is the highest and the conversion is limited. The chemical gradient is consequently low (68 wt. % to 60 wt. %). However, for y = 4 mm and higher, the fluid velocity decreases, the residence time and the heat flux increase (mass fraction difference at y = 4 mm of 48 wt. %). Thus, the pyrolysis is impacted by the border effect (wall of the cell). It is possible to plot it as a function of the x-axis, i.e. the porous medium thickness (Figure 13b). The Figure 13b shows that after 1 mm of porous flow, the wall impacts the flow on a distance of 1.5 mm. As a consequence, the radius of the permeable sample should not be lower than 1.5 mm. It should even be much more in order to limit its effect. The thicker the porous material in the present permeation test cell, the bigger should be the diameter in order to limit this border effect. This result is of major importance to design the experiments and to analyse the data. This is related to both dynamic and thermal boundary layers. Focusing on the porous zone, Figure 14 shows the profiles of pressure, of density, of temperature and of the axial velocity along the y-coordinate at the inlet and at the outlet surface of the porous medium. The pressure drop, due to the permeable sample, is 6800 Pa and it shows a uniform distribution on both sides of the porous zone. Although the pressure does not change in the crosswise direction (Figure 14a), the density decreases moving from the axis of symmetry toward the walls and from the upstream side to the downstream one (Figure 14b). This is directly due to the heat transfers (Figure 14c) and it acts on the fluid velocity (Figure 14d).

The spatial thermal gradient at the porous medium inlet presents a parabolic shape (Figure 14c) while it tends to be linear at the outlet. This could be due to the software assumptions considering equilibrium between solid and fluid. This may also be realistic in case of low speed flow (the conduction gets more important than the convection). In this latest case, it could be interesting to determine the length required by the fluid to present a complete linear profile. This would correspond to the thermal establishment length. Another interesting point is the location of the minimum velocity (Figure 14d). It is around y = 0.005 m at the porous medium inlet and y = 0.007 m at the outlet. This minimum is spatially shifted due to the thermal or chemical phenomena. Non-reactive hot flow could help understanding the reason of such a speed shift. Thus, a similar test case without reactions has been achieved for comparison purpose (Figure 15). Different temperature and velocity profiles have been obtained in comparison with Figure 14c,d. The difference is attributed to the chemistry, which impacts the molecular weight profile. The velocity profiles at the 0 mm and 3 mm positions are perfectly similar (the velocity amplitude only is shifted but not spatially). As a consequence, the modification of the minimum speed observed on

IV. Conclusion

The numerical simulation and the experimental tests are two complementary approaches, which are both suited for the study of fuel pyrolysis in porous media. Experimentally, all possible phenomena are taken into account with their intrinsic complexity and coupling, possibly under extreme conditions (1200 K and 6 MPa). Numerically, the heterogeneity profiles and the effect of the test bench itself (size and design of the permeation cell) have been observed finely despite the validity of the results is subject to the numerical, physical and chemical assumptions. A European Collaboration has been set up to provide a wide expertise in this field and to propose a validated numerical tool using detailed chemistry. A specific validation under no flow conditions has been first achieved for the chemistry since this is a major point in the numerical simulation of the reactive flow. After validation of the CFD simulation under inert conditions, its application to reactive test case presents average agreement with the experimental results. The computation method is judged promising. A modification by 25 % of the Forchheimer's permeability tends to confirm the accuracy which was suggested experimentally (±30 %). The pressure sensor position is suggested to be of importance considering the pressure drop measures. At least qualitatively, the relationships between the heat and mass transfers and the chemistry have been considered. The effects of the choice of the kinetic schemes on the multiphysics' simulation have been presented. Detailed mechanisms are clearly necessary. The comparison between numerical and experimental data appears difficult due to the complexity of accurately representing the experimental conditions. The numerical simulation showed strong benefits in understanding a possible border effect related to the permeable sample thickness. The relationship between the temperature field and the reaction zone tends to demonstrate the major impact of chemistry in the heat flux repartition. As a future work, other detailed kinetic mechanisms should be tested to reinforce the chemical modelling and to determine which level of accuracy is necessary for the chemistry to have a good representation of the phenomena. The final aim is to consider methane for rocket engines application.

2 C

 2 constant pressure (J.mol -1 .K -1 ) = inertial resistance factor (m -1 ) fluid / solid (kg.m -3 )

  numerical computations with three kinetic mechanisms at 900 K and 3.5 MPa after 0.7 s (a) and after 10 s (b) and at 1200 K and 3.5 MPa after 10 s (c) and dodecane pyrolysis as a function of time at 1200 K and 3.5 MPa (d).

Figure 3 .

 3 Figure 3. Schematic of the permeation test cell with porous medium and associated measures [25].

Figure 4

 4 Figure 4 Comparison of the numerical / experimental pressure drop as a function of the mass flow rate at 300 K and 2.5 MPa.

Figure 5

 5 Figure 5 Static pressure distribution with streamlines at 300 K and 2.5 MPa with nitrogen.

Figure 6

 6 Figure 6 Contour of axial velocity and streamlines for dodecane at 300 K and 3.5 MPa.

Figure 7 Figure 8

 78 Figure 7 Species mass fraction (wt. %) experimentally quantified upstream (HPC), in the porous medium (IPM) and downstream (LPC) during dodecane pyrolysis at 65 mg.s -1 , 3.5 MPa and 820 K.Table 2 Numerical test conditions for comparison with fluctuating experiments. Inlet fluid temperature (K) Inlet pipe temperature (K)

  12 H 26 mass fraction (a), temperature distribution (b), enthalpy of reaction (c) and molecular weight (d) for the two-steps mechanism at 1200 K, 3.5 MPa and 60 mg.s -1

Figure 11 .

 11 Figure 11. Species mass fraction and residence time along the axis of symmetry for the third case at 60 mg.s -1 (a) and the fourth case at 65 mg.s -1 (b).

Figure 12 .

 12 Dodecane mass fraction inside the system for case 1 (a) and case 4 (b) with residence time and mass fraction distribution along the symmetric x-axis (c).

  Porous medium thickness (mm) Influence distance of the wall (mm) b) Figure 13. a: C 12 H 26 mass fraction in the crosswise direction at the inlet and at the outlet of the porous medium at 1200 K, 3.5 MPa and 65 mg.s -1 and b: effect distance from the wall as a function of medium thickness.

  Figure 14d is due to the chemical reactions which modify the dynamic boundary layer. (a), density (b), temperature (c) and axial velocity (d) along y-direction at x = 0 mm and x=3 mm at 1200 K, 3.5 MPa and 65 mg.s -1 . (a) and axial velocity (b) along y-direction at x = 0 mm and x =3 mm without pyrolysis reactions in the porous medium (1200 K, 3.5 MPa, 65 mg.s -1 ).

  4 ,C 2 H 4 .

		80	Experiments [3]				80	Experiments [3]
	Mass Fraction (%)	20 30 40 50 60 70	This numerical work		Mass Fraction (%)	20 30 40 50 60 70	This numerical work
		10					10
		0					0
		H 2	C H 4 C 2 H 4 C 2 H 6 C 3 H 6 C 4 H 8 C 6 H 1 2 C 1 2 H 2 6	a)	H 2	C H 4 C 2 H 4 C 2 H 6 C 3 H 6 C 4 H 8 C 6 H 1 2 1 2 H 2 6 C	b)
				50
			Mass Fraction (%)	10 20 30 40
				0
				H 2	C H 4 C 2 H 4 C 2 H 6 C 3 H 6 C 4 H 8 C 6 H 1 2 C 1 2 H 2 6

Acknowledgments

This work was supported by the ESA-ESTEC, Contract no.: 3-12861/09/NL/PA. The authors would like to sincerely thank M. Bouchez, B. Lenaour and F. Travajno from MBDA-France, E. El Tabach and C. Strozzi from PRISME, D. Blanc and D. Courilleau from IUT Bourges, D. Davidenko from CNRS-ICARE Orléans, A. Ingenito from the University of Rome La Sapienza (Italy) and all the students involved in the projects head by I. Fedioun for their technical and scientific support involving the computations and the experiments.


J.-F Bercher

On escort distributions, q-gaussians and Fisher information

Keywords: nonextensive theory, escort distributions, Rényi-Tsallis entropy, Fisher information PACS: 02.50.Cw

Escort distributions are a simple one parameter deformation of an original distribution p. In Tsallis extended thermostatistics, the escort-averages, defined with respect to an escort distribution, have revealed useful in order to obtain analytical results and variational equations, with in particular the equilibrium distributions obtained as maxima of Rényi-Tsallis entropy subject to constraints in the form of a q-average. A central example is the q-gaussian, which is a generalization of the standard gaussian distribution.

In this contribution, we show that escort distributions emerge naturally as a maximum entropy trade-off between the distribution p(x) and the uniform distribution. This setting may typically describe a phase transition between two states. But escort distributions also appear in the fields of multifractal analysis, quantization and coding with interesting consequences. For the problem of coding, we recall a source coding theorem by Campbell relating a generalized measure of length to the Rényi-Tsallis entropy and exhibit the links with escort distributions together with pratical implications.

That q-gaussians arise from the maximization of Rényi-Tsallis entropy subject to a q-variance constraint is a known fact. We show here that the (squared) q-gaussian also appear as a minimum of Fisher information subject to the same q-variance constraint.

INTRODUCTION

Escort distributions have been introduced in statistical physics for the characterization of chaos and multifractals. These distributions P are a simple one parameter transformation of an original distribution p according to P(x) = p(x) q ´p(x) q dx .

The parameter q behaves as a microscope for exploring different regions of the measure p: for q > 1, the more singular regions are amplified, while for q < 1 the less singular regions are accentuated. This behavior is illustrated in Figure 1.

In Tsallis extended thermostatistics, the escort-averages, defined with respect to an escort distribution, have revealed useful in order to obtain analytical results and variational equations, with in particular the equilibrium distributions obtained as maxima of Rényi-Tsallis entropy subject to constraints in the form of a q-average. In applied fields, Tsallis distributions (q-distributions) have encountered a large success because of their remarkable agreement with experimental data, see [START_REF] Tsallis | Introduction to Nonextensive Statistical Mechanics[END_REF], [START_REF] Hilhorst | A note on q-Gaussians and non-Gaussians in statistical mechanics[END_REF], [START_REF] Vignat | Why is the detection of q-Gaussian behavior such a common occurrence?[END_REF], and references therein. Behavior of the mapping p(x) → P(x) = p(x) q / ´p(x) q dx. If U denotes the uniform distribution, it can be shown that I(U||P) ≥ I(U||p) if q > 1with the reverse inequality if , where I(p||q) denotes the Kullback-Leibler divergence from p to q. Examples of q-gaussian distributions with γ = 1 and σ 2 = 1, for several values of q. In the first case, q ≤ 1 , the distributions have compact support. In the second case where q ≥ 1, the distributions have an infinite support and heavy tails. The limit case q = 1 is the standard gaussian distribution.

A central example is the q-gaussian distribution, which has the following probability density function

f (x) = 1 Z q (β ) 1 -(1 -q)β x 2 γ 1-q + , q = 1 (2)
where γ and β are positive parameters, and where we use the notation (x) + = max {x, 0} . The partition function converges for q < 2γ + 1, and the variance exists for q < 2γ/3 + 1.

Examples of q-gaussian distributions for different values of q are given on Figure 2, with γ = 1 and with variance σ 2 = 1.

In this contribution, we show that escort distributions (1) emerge naturally as a maximum entropy trade-off between the distribution p(x) and the uniform distribution. This setting may typically describe a phase transition between two states. Then, looking for the distribution p with maximum entropy subject to a constraint on the mean energy computed with respect to the escort distribution, we arrive at the q-gaussian (2).

But escort distributions also appear in other fields with interesting consequences. In the problem of coding, we recall a source coding theorem by Campbell relating a generalized measure of length to the Rényi-Tsallis entropy. We show that the associated optimal codes can easily be obtained using considerations on escort-distributions and that this provide an easy implementation procedure. We also show that these generalized lengths are bounded below by the Rényi entropy.

That q-gaussians arise from the maximization of Rényi-Tsallis entropy subject to a qvariance constraint is a known fact. We show here that the q-gaussian (with γ = 2) also appears as a minimum of Fisher information subject to the same q-variance constraint, and of course, that we recover the standard results in the limit q = 1.

ESCORT DISTRIBUTIONS Escort as a model of displaced equilibrium

It has been noted that Tsallis' extended thermodynamic seems appropriate in the case of modified, perturbated, or displaced classical Boltzmann-Gibbs equilibrium. This means that the original MaxEnt formulation "find the closest distribution to a reference under a mean constraint" may be amended by introducing a new constraint that displaces the equilibrium, as we discussed in a previous MaxEnt workshop [START_REF] Bercher | An amended MaxEnt formulation for deriving Tsallis factors, and associated issues[END_REF]. The partial or displaced equilibrium may be imagined as an equilibrium characterized by two references, say r and q. Instead of selecting the nearest distribution to a reference under a mean constraint, we may look for a distribution p * simultaneously close to two distinct references: such a distribution will be localized somewhere 'between' the two references r and q. For instance, we may consider a global system composed of two subsystems characterized by two prior reference distributions. The global equilibrium is attained for some intermediate distribution, and the observable may be, depending on the viewpoint or on the experiment, either the mean under the distribution of the global system or under the distribution of one subsystem. This can model a fragmentation process: a system Σ(A, B) fragments into A, with distribution f , and B with distribution g, and the whole system is viewed with distribution p * that is some intermediate between f and g. This can also model a phase transition: a system leaves a state g toward f and presents an intermediate distribution p * . This intermediate distribution shall minimize its divergence to the initial distribution g(x), but also be 'not too far' from the attractor distribution f (x). This can be stated as I(p||g) -I(p|| f ) ≤ θ , or equivalently as I(p|| f ) ≤ θ ′ , where I(p||q) denotes the Kullback-Leibler divergence from a distribution p to another distribution q. Remark that the first constraint can be interpreted as a constraint on the mean log-likelihood. The problem can be written as follows:

   min p I(p||g) s.t. I(p||g) -I(p|| f ) ≤ θ or I(p|| f ) < θ ′ , (3)
where "s.t." stands for "subject to". The optimum distribution solution of these two equivalent problems is

f * (x) = f (x) q g(x) 1-q ´f (x) q g(x) 1-q dx , (4)
The solution of the minimization problem satisfies a Pythagorean equality: I(p||g) = I(p|| f * ) + I( f * ||g) for any distribution p such that θ = I(p||g) -I(p|| f ). It is interesting to note that the solution ( 4) is nothing but a generalized version of the escort or zooming distribution of nonextensive thermostatistics. Obviously, one recover (1) when g(x) is the (improper) uniform measure. With the expression of the solution p * (x), we obtain that

I( f * ||g) = qθ -log( ˆf (x) q g(x) 1-q dx) = qθ -(q -1)I q ( f ||g), (5)
where we have recognized the Rényi divergence with index q, I q ( f ||g) from f to g. Suppose now that we get an additional constraint, in the form of an observable, given as some expectation under p * , the equilibrium distribution between f and g. Then, a new question arises, that is the determination of the more general distribution f , compatible with this constraint.

We may use a MaxEnt point of view that consists in further minimizing the Kullback-Leibler information divergence I( f * ||g), with respect to f . This finally amounts to the minimization of the Rényi divergence in (5), I q ( f ||g) subject to the mean constraint. In the case of a uniform distribution for g and a variance constraint we obtain the so-called q-variance constraint

σ 2 q = ´x2 f (x) q dx ´f (x) q dx . (6)
Obsiously, the problem reduces to the maximization of Rényi-Tsallis entropy subject to this q-variance constraint, whose solution is known to be the standard q-gaussian distribution (2) with γ = 1.

Another possibility is to look for the less informative distribution in the Fisher sense, subject to the same constraint. We will see that this also leads to a q-gaussian distribution.

Rényi entropy, escort distributions and multifractals

Tsallis was originally inspired by the multifractals for the proposal of a generalized entropy. In this setting, one describes the properties of a singular measure in terms of singularity spectrum and generalized dimensions. A standard construction is to imagines that the support of the measure is covered with boxes of size l. Then if p i is the probability in the i-th box, one puts p i (l) ∼ l -α i , where α i is the singularity exponent. Similarly, the number of boxes N(α) with an exponent α defines the singularity spectrum f (α) with N(α) ∼ l -f (α) . An alternative description is provided by the q-th 'moments' of the measure N q = ∑ i p i (l) q ∼ l (q-1)D q ∼ l τ(q) , where D q is the 'Rényi dimension' and τ(q) is the correlation exponent. The important connection between the two descriptions is that the pairs (α, f (α))and (q, τ(q)) are related by a Legendre transform [START_REF] Riedi | Multifractal processes[END_REF], i.e. f (α) =qα -τ(q) with α = dτ(q)/dq.

It is well known that MaxEnt problems, as well as the standard thermodynamics are characterized by a beautiful Legendre structure. Thus, we readily see that the previous Legendre transform can be given a maximum entropy interpretation, where f (α), the singularity spectrum, is an entropy and τ(q) is the free energy. Here, the free energy τ(q) has the form τ(q) = log ∑ i p i (l) q = log ∑ i exp(q log p i (l)), and we may identify τ(q) as τ(q) = log ∑ i exp(β E i ) with β = q the inverse temperature and E i = log p i (l) the energy (also called here the 'bit-number'). This leads us to an easy thermodynamic interpretation: exp(β E i ) is the Boltzmann factor which results from the maximization of the entropy -∑ i P i log P i subject to the mean constraint

U = α = E P [E i ] = E P [log p i (l)].
The standard resolution of this problem [START_REF] Harte | Multifractals: Theory and Applications[END_REF][START_REF] Bercher | Tsallis distribution as a standard maximum entropy solution with 'tail' constraint[END_REF] leads us to the solution

P i = p q i ∑ i p q i
, which is nothing else but the escort distribution of order q of {p i }. Obviously, we then obtain that the singularity spectrum is the entropy of the escort, and the mean exponent is the mean of the bit-number:

f (α) = ∑ i P i log P i and α = E P [log p i (l)].
These relationships are the rationale basis for the method introduced in [START_REF] Chhabra | Direct determination of the f (α) singularity spectrum[END_REF] for the direct estimation of the singularity spectrum. In the latter paper, the escort distributions were introduced as a simple operational tool. They were later named "escort distributions" in the book by Beck [START_REF] Beck | Thermodynamics of Chaotic Systems[END_REF].

Source coding

We also have an interesting connection with source coding, which was first reported in [START_REF] Bercher | Source coding with escort distributions and Rényi entropy bounds[END_REF]. The mains formulas and comparisons are given on Table 1.

In source coding, one considers a source that produces a set of symbols X = {x 1 , x 2 , . . . x N } with probabilities p i . The aim of source coding is to map each symbol x i to a codeword c i of length l i expressed using the D letters of the alphabet. The Shannon source coding theorem (noiseless coding theorem) indicates that the expected length L of any uniquely decodable code is bounded below by the entropy of the source, H 1 (p).

It is well-known that Huffman coding yields a prefix code which minimizes the expected length and approaches the optimum limit l i =log D p i for the individual symbols. Other forms of lengths have also been considered, the first and fundamental contribution being Campbell's one [START_REF] Campbell | A coding theorem and Rényi's entropy[END_REF]. In Shannon's result, low probabilities yield very long words. In order to penalize the apparition of such long codes, Campbell proposed a different length measure, featuring an exponential account of the elementary lengths of the codewords.

The remarkable result [START_REF] Campbell | A coding theorem and Rényi's entropy[END_REF] is that just as Shannon entropy is the lower bound on the average codeword length of an uniquely decodable code, the Rényi entropy of order q, with q = 1/(β + 1), is the lower bound on the exponentially weighted codeword length, see line 2 of Table 1.

In order to increase the impact of the longer lengths with low probabilities, the Campbell's length uses an exponential of the length. A different approach to the problem can be to modify the weigths in the linear combination, so as to raise the importance of the terms with low probabilities. A simple way to achieve this is to deform, flatten, the original probability distribution so as to scale down its features with highest probabilities and raise the areas of lower probability. Of course, a very good candidate is the escort distribution, which leads us to the 'average length measure' M q = ∑ N i=1 P i l i . For the virtual source with distribution P, the standard expected length is M q , and the classical Table 1. Source coding with generalized lengths, with P i = p q i / ∑ i p q i and where log D denotes the base D logarithm length bound opt. length

Shannon L = ∑ i p i l i H 1 (p) l i = -log D (p i ) Campbell C β = 1 β log D ∑ N i=1 p i D β l i H 1/(β +1) (p) l i = -log D (P i ) Gen. mean M q = ∑ N i=1 P i l i H 1 (P) l i = -log D (P i ) Mixed L q = 1 q-1 log D ∑ N i=1 p q i ∑ j p q j D (q-1)l i H q (p) l i = -log D (p i )
Shannon noiseless source coding theorem immediately applies, leading to M q ≥ H 1 (P) with equality for l i =log D P i , which is nothing else but Campbell's optimal lengths. The simple relation l i =log D P i for the minimization of M q has an immediate application for the computation of Campbell's lengths: it suffices to feed a standard coding algorithm, namely a Huffman coder, with the escort distribution P instead of the natural distribution p. Another measure of length, Table 1, line 4, can mix both a an exponential weight of individual lengths l i and an escort distribution. This length of order q is lower bounded by the Rényi entropy of the same order. Interestingly, the optimum legths are nothing but the optimal lengths in the Shannon coding theorem.

FISHER, ESCORT AND q-GAUSSIANS

Let f denote the probability density of a random variable X. Then the Fisher information with respect to a translation parameter is defined as

I[ f ] = ˆ d ln f (x) dx 2 f (x)dx = ˆ d f (x) dx 2 1 f (x) dx, (7)
where f (x) is supposed differentiable and both f (x) and its derivative f ′ (x) are square integrable on R.

The importance of the Fisher information as a measure of the information about a parameter in a distribution is well known. It has intricate relationships with maximum likelihood and has many implications in estimation theory, as exemplified by the Cramér-Rao bound which provides a fundamental lower bound on the variance of an estimator. It is used as a method of inference and understanding in statistical physics and biology, as promoted by Frieden [START_REF] Frieden | Science from Fisher Information: A Unification[END_REF]. Information theoretic inequalities involving Fisher information have attracted lot of attention for characterizing statistical systems through their localization in information planes, e.g. the Fisher-Shannon information plane [START_REF] Vignat | Analysis of signals in the Fisher-Shannon information plane[END_REF][START_REF] Romera | The Fisher-Shannon information plane, an electron correlation tool[END_REF] or the Cramér-Rao information plane.

If we return to the problem of selecting the distribution f (x) with a given q-variance, as indicated in subsection , we see that such distribution can be sought in a minimum Fisher sense. Actually, there are at least two possibilities that makes sense: (i) find f (x) with a given variance such that its escort has minimum Fisher information, or (ii) find a minimum Fisher information f (x) such that its escort has a given variance, i.e. f (x) has a given q-variance. Due to lack of space, we will only consider here the second variational problem:

     inf f I [ f ] , s.t. σ 2 q = ´x2 f (x) q dx ´f (x) q dx , s.t. f (x) ≥ 0, ´f (x)dx = 1. (8)
The solution to (i) can be simply deduced from the solution to (ii).

The main difficulty is that even though the Fisher information I[ f ] is a strictly convex function, the set defined by the q-variance constraint is not convex and uniqueness of a solution is not guaranteed. It is convenient to introduce a new parameter and reformulate the awkward q-variance constraint into two simpler elementary constraints, that is N q = ´f (x) q dx and V q = ´x2 f (x) q dx, with, of course, σ 2 q = V q /N q . Then, we can note that the initial problem is equivalent to the following two steps procedure:

         inf N q      inf f I [ f ] , s.t. N q = ´f (x) q dx and V q = σ 2 q N q = ´x2 f (x) q dx s.t. ´fN q (x)dx = 1 . ( 9
)
provided that the second problem is feasible with the parametric solution obtained in the first step. For now, we only consider the first step. Again and as in [START_REF] Chhabra | Direct determination of the f (α) singularity spectrum[END_REF], although the objective functional I[ f ] is strictly convex, since the constraints sets are not convex then the uniqueness of a minimum cannot be guaranteed. However, it is still possible to relax the equality constraints into inequalities, e.g. Sign(q -1) ´f (x) q dx -N q ≤ 0 so that the problem becomes a convex optimization problem, whose solution is known to be unique. If the constraints are active at the solution, that is with equality, then we obtain the solution of the initial problem using the (relaxed) convex formulation. Let a and b denote the Lagrange multipliers associated with the inequality constraints. These multipliers are necessarily nonnegative. Now, can check that with the following candidate f

(x) = u(x) 2 = α 2 1-q 1 -(1 -q)β x 2 2 1-q (10)
satisfies the Euler-Lagrange equation associated with the relaxed problem, with a = 4Sign(q -1) α 2 β q and b = -4Sign(q -1) α 2 β 2 (1-q) q . With β > 0, we check that these Lagrange multipliers are simultaneously positive for q > 1, but do not have the same sign for q < 1. This means that for q > 1, the q-gaussian (10) is a possible solution of the relaxed problem.

With a, b positive the Lagrangian is a convex function of f and of f ′ (I[ f ] is a strictly convex function of f and that this is also true pointwise). In that case, any singular point corresponds to the global minimum [15, theorem 2.1]. Therefore, we obtain that the q-gaussian solves the relaxed minimization problem. Finally, it always possible to choose α so that the solution is normalized, then select β and adjust the parameter N q so that the inequality constraints are active at the solution, i.e. are verified with equality. As a result, we finally obtain the global minimum f * of the initial problem is obtained with the q-gaussian (2) with γ = 2. This approach does not cover the case q < 1. In this case, we also obtain that a q-gaussian, with compact support, is the minimizer of the Fisher information, but the proof of the general case is more involved and omitted here. A direct consequence of this result is a generalization of Cramér-Rao inequality which readily follows from the inequality I[ f ] ≥ I[ f * ] for all f with a given q-variance.

CONCLUSIONS AND FURTHER WORK

In this contribution, we have underlined the apparition of escort distributions in several domains. We have proposed a view of escort distributions as a trade-off between two distributions, in the context of displaced equilibriums or phase transition. When q varies, the escort distribution defines a path between two distributions, that shall be characterized; e.g. the thermodynamic length which is defined as the flux of Fisher information along the path. Associated with escort distributions, we recall that the q-gaussians are Rény-Tsallis Maxent distributions. This characterization is completed by the fact that q-gaussians are also minimum Fisher information distributions. Further work includes the study of the interaction between (generalized) Fisher information and Rényi-Tsallis entropy, as well as the study of the Fisher path, defined as the curve of minimum Fisher distributions when q varies.

Figure 1 .

 1 Figure1. Behavior of the mapping p(x) → P(x) = p(x) q / ´p(x) q dx. If U denotes the uniform distribution, it can be shown that I(U||P) ≥ I(U||p) if q > 1with the reverse inequality if , where I(p||q) denotes the Kullback-Leibler divergence from p to q.

Figure 2 .

 2 Figure2. Examples of q-gaussian distributions with γ = 1 and σ 2 = 1, for several values of q. In the first case, q ≤ 1 , the distributions have compact support. In the second case where q ≥ 1, the distributions have an infinite support and heavy tails. The limit case q = 1 is the standard gaussian distribution.


Laurent Fiack

Benoît Miramond

Nicolas Cuperlier

FPGA-based vision perception architecture for robotic missions

Most of the robot behaviors are based on visual sensing and perception. This paper describes a smart camera composed of a full-hardware vision architecture coupled with an embedded camera sensor. The hardware architecture corresponds to low-level visual perception processes. The integration of such a system onto the robot enables not only to accelerate the visual processing till real-time behavior but also to compress the data-flow at the output of the camera. The results obtained during indoor robotic missions show an important reduction factor of data communication.

I. INTRODUCTION

Real-time visual processing represents a major challenge for autonomous robots. Most of the robot behaviors are based on visual sensing and perception: navigation, object recognition and manipulation, target tracking, and even social interactions with humans. The visual systems require large computing capabilities which make them hard to embed. Smart cameras are proposed as an alternative to get visual low-level processing back into the robot.

These visual systems are not considered isolated anymore but as part of an architecture integrated in its environment. They take into account several parameters related to the dynamic properties of the systems they belong to (see active vision [START_REF] Ballard | Animate vision[END_REF]). These visual processing algorithms strongly depend on the dynamics of interactions between the system and its environment by continuous feedbacks regulating even the low-level visual stages such as the attentional mechanisms in biological systems.

The robotic missions considered in this paper consists in a subset of a complex cognitive system allowing a robot equipped with a charge-coupled device (CCD) camera to navigate and to perceive objects. The global architecture in which the visual system is integrated is biologically inspired and based on the interactions between the processing of the visual flow and the robot movements (Per-Ac architecture [START_REF] Gaussier | Perac: a neural architecture to control artificial animals[END_REF]). The learning of the sensorimotor associations allows the system to regulate its dynamics [START_REF] Maillard | Perception as a dynamical sensori-motor attraction basin[END_REF] and, therefore, navigate, recognize objects, or create a visual compass.

In this paper, we aim at designing an embedded visual processing system in the form of a single chip that could be used for building the CCD-based smart camera of our robot. On one hand, the embedded processing part should be configurable in order to allow a variety of navigation missions. On the other hand, the architecture should also provide intensive computation capabilities to deal with lowlevel image processing. Finally, the system should reduce the amount of data communication at the output of the camera without loss of information.

Thus, the proposed architecture combines video sensing, perception processing and communication into a System-on-Chip (SoC) embedded onto the robot. This vision system is designed as a full hardware architecture deployed onto a FPGA device. The system is able to process video frames in real-time and communicates intensively with an embedded computer running high-level cognitive tasks that depend on the robot mission. GPU-based solutions provide similar performances but at the cost of important energy consumption that make them unusable in an embedded context.

The rest of the paper is organized as follows. Section II presents the context of this work, especially it describes the different types of robotic missions. For each mission, the environment and the cognitive tasks are different. For that reason the robot needs specific visual features that are provided by the configurable vision SoC. Section III describes the architecture of this system and the parameters that are used to configure the smart camera. Section IV presents the experimental platform and provides the results obtained during robotic missions. Finally, we conclude and outline future works in section V.

II. CONTEXT

Among different perceptual modalities, vision is certainly the most important building block of a bio-inspired cognitive robot. Vision enables robots to perceive the external world in order to perform a large range of tasks such as navigation, object tracking and manipulation, and even interaction with humans. This broad range of tasks often relies on different models and architectures. According to a given task, the visual processing often used only a subset of the available algorithms like optical flow, feature points extraction and recognition (over one or several scale spaces). For instance, navigation often relies on low resolution images to only capture the main regularities from the environment, whereas object recognition may need to characterize an object over several scales in order to take into account more details of the object.

The current visual system integrates a multiscale approach to extract the visual primitives. In doing so, it also allows a wider range of applications. Roughly the visual system provides a local characterisation of the keypoints detected on the image flow of an 8-bit gray-scale CCD camera. This local characterisation feeds a neural network which can associate motor actions with visual information: this neural network can learn, for example, the direction of a displacement of the robot as a function of the scene recognition. The studied visual system can be divided into two main modules:

• a multiscale mechanism for characteristic points extraction (keypoints detection), • a mechanism supplying a local feature of each keypoint.

A. The multiscale keypoints detection

The multiresolution approach is now well known in the vision community. A wide variety of keypoints detectors based on multiresolution mechanisms can be found in the literature. Amongst them are the Lindeberg interest point detector [START_REF] Lindeberg | Feature detection with automatic scale selection[END_REF], the Lowe detector [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] based on local maxima of the image filtered by difference of Gaussians (DoGs) or the Mikolajczyk detector [START_REF] Mikolajczyk | Scale & affine invariant interest point detectors[END_REF], where keypoints correspond to those provided by the computation of a 2D Harris function and fit local maxima of the Laplacian over scales. The visual system described here is inspired from cognitive psychology. The used detector extracts points in the neighbourhood of the keypoints, which are sharp corners of the robots visual environment. More precisely, the keypoints correspond to the local maxima of the gradient magnitude image filtered by DoGs. The detector is characterised by a good stability. Following the detector itself, the system provides a list of sorted local features by the way of competition between the corresponding keypoints.

Keypoints are detected in a sampled scale space based on an image pyramid. Pyramids are used in multiresolution methods to avoid expensive computations due to filtering operations. The algorithm used to construct the pyramid is detailed and evaluated in [START_REF] Crowley | Fast computation of scale normalised gaussian receptive fields[END_REF]. The pyramid is based on successive image filtering with 2D Gaussian kernels normalised by a factor S.

These operations achieve successive smoothing of the images. Two successive smoothing are carried out by two Gaussian kernels with variance σ 2 = 1 and σ 2 = 2. The scale factor doubles (achievement of an octave) and thus the image is decimated by a factor of two without loss of information. The same Gaussian kernels can be reused to continue the pyramid construction. Interestingly, the kernel sizes remain small (half-width and half-height of 3σ) allowing a fast computation of the pyramid. Finally, the images filtered by DoGs in the pyramid can be simply obtained by subtracting two consecutive images. Keypoints detected on the images are the first N local maxima existing in each DoG image of the pyramid. Thus, the keypoint research algorithm orders the N first local maxima according to their intensities and extracts their coordinates. The shape of the neighbourhood for the research of maxima is a circular region with a radius of 21 pixels. This value is one of the system parameters. The number N which parametrises the algorithm corresponds to a maximal number of detections. Indeed, the robot may explore various visual environments (indoor versus outdoor) and particularly more or less cluttered scenes may be captured (e.g., walls with no salience versus complex objects). A detection threshold (γ) is set to avoid nonsalient keypoints. This threshold is based on a minimal value of the local maxima detected. The presence of this threshold is even more important in the lowest resolutions since the information is very coarse at these resolutions. This particularity of the algorithm confers it a dynamical aspect. Precisely, the number of keypoints (and consequently the number of local features) depends on the visual scene and is not known a priori. Furthermore, the threshold γ could be set dynamically through a context recognition feedback but discussing here this mechanism is not our current purpose. However, even if this threshold is considered as a constant value, the number of detected keypoints varies dynamically according to the input visual scene. Consequently, the number of computations (neighbourhood extractions) also depends on the input data.

B. The local image feature extraction

At this stage, the neighbourhood of each keypoint has to be characterised in order to be learnt by the neural network. Existing approaches to locally characterise keypoints are numerous in the literature: local jets, scale invariant feature transform (SIFT) and its variants, steerable filters, and so forth. In the current application, we simply reuse a viewbased characterisation where keypoint neighbourhoods are represented in a log-polar space. This representation has good properties in terms of scale and rotation robustness. The local feature of each keypoint is, therefore, a small image resulting from the log-polar transformation of the neighbourhood. The neighbourhoods are extracted from the gradient magnitude image at the scale the keypoint was found by the detector. Each neighbourhood extracted is a ring of radius (5, 16) pixels. Excluding the small interior disc avoids multiple representations of the central pixels in the log-polar coordinates. The angular and logarithmic radius scale of the log-polar mapping are sampled with 16 values. Each feature is thus an image of dimension 16x16 pixels. The sizes of the rings and feature images have been determined experimentally for an indoor object recognition. The given parameters represent a tradeoff between stability and specificity of the features. Finally, the small log-polar images are normalised before their use by the rest of the neural architecture. By associating the data provided by the visual system with actions, the global system allows the robot to behave coherently in its environment [START_REF] Maillard | Perception as a dynamical sensori-motor attraction basin[END_REF]. Generally, the visual system must not be considered isolated but integrated in a whole architecture whose modules interact dynamically with each other and through the environment. Hence, the evaluation of the parameters of the visual system depends on the rest of the robot system architecture.

III. HARDWARE ARCHITECTURE

This section describes the hardware implementation of the bio-inspired visual system, which feeds the neural network of the mobile robot. The visual system is implemented into a FPGA device which communicates through Ethernet with an embedded computer running the neural network (Figure 1).

This vision system requires heavy computing power and can't be embedded as software respecting real-time constraints (the frame-rate of the camera). That's why we realized a hardware implementation through an algorithm/architecture adequacy process [START_REF] Verdier | Using high-level rtos models for hw/sw embedded architecture exploration : Case study on mobile robotic vision[END_REF]. Integrating these first steps of vision perception into the robot itself brings interesting reductions of the data communication. The camera does not act as a primitive sensor anymore but as a smart camera modeling pre-attentional mechanisms that only keeps the most salient visual points. The Gaussian pyramid is presented in figure 2.

A. Architecture of the pyramid

The gradient and the Gaussian filters are the results of convolutions by a kernel. As often in hardware vision, the pixels of an image come as a flow, and they are processed with the same rhythm. Instead of moving the kernel over the image, the image moves over an architecture which computes the matrix product. A number of pixels must then be stored in first-in first-out (FIFO) memories, depending on the kernel size and the image width. The output pixels are produced with a fixed latency, depending also on these two parameters.

The DoG just computes the difference of it's two inputs. Because of the latency of the Gaussian filters, a FIFO is needed for one of its inputs. The subsampling decimates the input image by a factor of two, and thus divides the coordinates of the pixels by two.

In the first part of his project, Lefebvre [START_REF] Lefebvre | Phd report : Exploration architecturale pour la conception d'un systeme sur puce de vision robotique, adequation algorithmearchitecture d'un systeme embarque temps reel[END_REF] proposed an interface for its processing IPs. With this interface, the input and output buses are composed of the pixels value and their coordinates. A standard interface allows for modularity of these IPs. It is then easy to assemble the processing IPs (gradient, gauss, dog and subsampling) to build the hardware Gaussian pyramid (Figure 2).

B. Keypoint search and sorting

The next step consists in searching the local maxima of a DoG output, to sort them by magnitude, and finally to keep the N most significant points. Each pixel is compared to all the pixels in the neighborhood disc. To follow the pixel rhythm, this comparison must be done in one clock cycle. Once a potential keypoint is found, the exclusion disc avoids two keypoints to be to close.

These points then enter in the sort IP. This component consists in a register bank which stores the value, the coordinates, and an index corresponding to the IP that will process the log-polar mapping (Figure 4). These points are stored at a position depending on its value. When a point must be inserted between two registers, the lower part can be shifted down.

C. Neighborhood extraction and log-polar mapping

Two IPs are responsible for the log-polar mapping (Figure 4), Address Generator and Transform. The Address Generator converts the Cartesian address into the log-polar one. A test ensures that the incoming flow is in the neighborhood disc of the keypoint, and then its coordinates are removed from the coordinates of the flow. These coordinates are the input of a look-up table (LUT) which stores the final logpolar coordinates.

The Transform IP reads the incoming pixels, and stores them into a memory at the address computed by the Address Generator. This IP is built as a double-buffer memory to compute an image by reading the previous one. The zones of the outer rings containing more than one pixels must be averaged. An accumulator is then added to the memory. As the address conversion is done by a LUT, the values for the average division can also be stored into a LUT.

We need to duplicate these two IPs as many as the number of keypoints we want to detect (Figure 4). The keypoints are dispatched to an Address Generator depending on the index coming from the sort IP. As we cannot shift down the content of the entire memory of a Transform bloc, we can only overwrite the memory corresponding to the lower keypoint. That's why we need a connection between the rank of the keypoint and its index. Finally, a CPU can read all the data by selecting an index. A memory (not represented) stores the value of the keypoint, its coordinates, and the relation between its rank and its index, so that the CPU doesn't have to sort the values again. It can read the log-polar feature as a specific RAM memory.

D. Connection with the PC

We choose to use Ethernet and TCP/IP as it offers sufficient throughput, and was already working on the neural network software hosted by the PC. We use a NIOS II soft-processor from Altera to manage the TCP/IP stack. The Ethernet MAC is built has a hardware peripheral. We designed a custom memory-mapped peripheral to read the keypoints values, their coordinates and their associated logpolar features. The global architecture of the smart-camera is depicted in figure 5 .

5. SOPC. The logpol reader makes the junction between the visual system and the microprocessor world.

IV. EXPERIMENTS AND RESULTS

A. Vision quality and architecture sizing

The first stage of the pyramid have been tested, and works in real-time. The architecture supports the flow of 320 × 240 images at 22 frames-per-second given by the image sensor.

The synthesis results are given in Figure 6 for the Altera Cyclone IV EP4CE115 FPGA device of the Terasic DE2-115 board. The FPGA architecture is composed of 114480 LEs (Logic Elements) and 3981312 Memory bits. These results are given for the first stage of the pyramid, for 320 × 240 images, with a 16-bits bus width. The radius for the research, the exclusion and the extraction is set to 21 pixels, and the radius of the blind zone is set to 5 pixels. The architecture keeps a maximum of N = 16 keypoints. The Capture represents all the hardware responsible for reading the images. The DoG represents the two Gaussian filters, and the DoG itself. The Logpol is the IP shown in Figure 4. The SOPC contains the NIOS II, the Ethernet MAC, the SDRAM controller and the specific peripheral which reads the values from Logpol.

IP

The execution of a single stage of the pyramid already uses almost the full FPGA resources. Working with the entire pyramid would need more recent and denser FPGA matrices.

B. Performances and communication

The architecture generates 16 features of 16 × 16 pixels of 16 bits for each 320 × 240, 8-bits frame. We can easily compute the compression ratio : This architecture allows a much smaller computing time than a classical computer can offer for less power consumption. The visual system itself reduces the amount of data to transmit. The system is currently embedded onto the robot as depicted in figure 7. Moreover we are now working at optimizing the system dimensioning, especially the size of the Logpol IP. The size of the Search can also be reduced since it decreases quadratically with the research radius. But globally there will always be a trade off between the number of scale to compute and the resolution of the features to extract.

V. CONCLUSIONS

Visual processing in the context of autonomous robots is often a challenge both for real-time interactions with the environment and for communication bottlenecks with distant robotic platforms. This paper describes a FPGA-based smart camera executing low-level visual perception tasks directly onto the robot. The system composed of a CCD camera and an FPGA device only outputs a set of keypoints at each frame sampled by the camera. The keypoints detector is implemented as a full-hardware architecture. The system works in real-time at 22 frames per second and enables a compression ratio above 9 compared to a classical camera.

The visual system integrates a multiscale approach to extract the visual primitives. In this paper, we provided results for the first stage of the pyramid (high frequencies). We are now working at integrating the full system into more recent FPGA families with denser logic matrix.

Fig. 1 .

 1 Fig. 1. Schematic view of the robot

Fig. 2 .

 2 Fig. 2. Gaussian Pyramid. Each output corresponds to a frequency band of the input image.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Log-polar mapping. The image on the left represents the log-polar mapping with ρ = 5 and θ = 16, the image on the right represents the sampled mapping over the cartesian grid with an extraction radius of 20 and a blind zone radius of 4 pixels, and the third image is the result of the remapping in the (ρ, θ) space.

Fig. 7 .

 7 Fig. 7. The smart camera embedded onto the robot