Massyl Lagha
email: massyl.lagha@isae.fr

Sébastien Duplaa

Nicolas Binder

Xavier Carbonneau

Unified classification and characterization of axial turbomachines and propellers

Keywords: turbomachinery, characteristic map, (ψ, φ) formalism, performance analysis, fan, propeller

Reaching the mid/long-term air transport emission reduction goals imposed by both European and American standards impose increasing the propulsive systems' adaptability to various operating conditions, in order to maximize the aircraft overall efficiency all along the flight mission. This implies the enlargement of the design space of propulsive systems such that it can even be operated equally as a compressor or turbine, which leads to rethink the paradigm of designing turbomachines. The continuity in the definition and characterization of different types of turbomachines should be restored which is proposed through this contribution. Analytical relationships allowing to switch between compressor map, propeller map and (ψ, φ) map are developed. To minimize the inputs of the maps' conversion relations, a methodology to extract mean flow features from any characterization map is presented, namely the rotor outlet relative flow angle and mean streamline radius. The application of the characterization maps' conversion relations on a turbofan's single-stage axial fan and on a propeller allowed their validation through the physical coherence of the results. The flow features extraction methodology also showed very satisfying results with comparison to experiments. Eventually, the ability of the (ψ, φ) formalism as a powerful performance analysis tool for all kind of turbomachines is stressed out, which makes it the best candidate for the unified treatment of turbomachines.

Introduction

The challenging air transport emission reduction goals imposed by both the European (ACARE) [START_REF]Flightpath 2050 Europe's Vision for Aviation[END_REF] and the American (NASA) [2] standards require development of innovative propulsion technologies. In addition to the increase of efficiency that has to be achieved on propulsive systems, one must consider global optimization of flight operations. This optimization implies that the operating conditions of the propulsive system may have to vary within a wide range during the flight mission. The recent research tendencies head toward the introduction of geometrical variabilities on current propulsive systems [START_REF] Zarati | Engine performance and surge margins optimization by means of nozzles variability[END_REF][START_REF] Joksimović | Local and global analysis of a variable pitch fan turbofan engine[END_REF]. By 2050, the IATA [START_REF]IATA Technology Roadmap[END_REF] demonstrate that emission reduction goals cannot be achieved without shifting research towards disruptive technologies such as distributed propulsion.

The purpose of introducing geometrical variabilities or distributing propulsion is to propose a larger design space in order to increase the propulsive systems' adaptability to various operating conditions. This would allow considering maximization of aircraft overall efficiency all along the flight mission. In case of distributed propulsion, it would make it possible to partially recover power through the windmilling of a limited number of propulsive modules during descent phase. In order to enlarge design space and to propose dual operating design (compressor-turbine), one must reconsider the multiple turbomachines' classifications and go back to their most general definition, namely: "a mechanical device exchanging work with a continuous fluid flow through one or multiple rotors". This vision is necessary as the turbomachines' community has been divided since its early birth due to the intimate link between the field of turbomachinery and engineering companies. Theses companies make designs more specific and less public, which prevent development of a unified vision and tools for turbomachines' analysis. One of the most striking examples that illustrates this division is the very distinct treatment that is made for propellers that are considered as a specialized field of external flows aerodynamics while all other turbomachines are related to internal flows aerodynamics. For this reason, propellers are often not even considered as turbomachines whereas they perfectly match the abovementioned definition. The main differences between propellers and the other turbomachines lies in the lower solidity and work input of propellers which allow their treatment through more classical fluid mechanics tools that are usually attributed to external flows aerodynamics [START_REF] Wislicenus | Fluid Mechanics of Turbomachinery[END_REF]. Still, it is not a sufficient reason to break continuity between propellers and the rest of turbomachines.

In fact, as it can be depicted in many excellent reference literature [START_REF] Lakshminarayana | Fluid dynamics and heat transfer of turbomachinery[END_REF][START_REF] Dixon | Fluid mechanics and thermodynamics of turbomachinery[END_REF], turbomachines are firstly classified into different types and then treated independently using different design and analysis tools. Each type is also restricted to a precise specification window as presented by Japikse and Baines [START_REF] Japikse | Introduction to turbomachinery[END_REF] (cf. figure 1).

This compartmentalization is mainly due to the work of Cordier [START_REF] Cordier | Ähnlichkeitsbedingungen für Strömungsmaschinen[END_REF] further developed by Balje [START_REF] Balje | Turbomachines : a guide to design selection and theory[END_REF] aiming to correlate non-dimensional numbers (specific speed ω s and specific diameter D s ) to the optimal turbomachine's geometry. Those non-dimensional numbers, as defined by Balje are:

ω s = ω ( V ) 1 /2 (∆h i ) 3 /4 (1)
D s = D (∆h i ) 1 /4 ( V ) 1 /2 (2)
Cordier and Balje computed these numbers from measurements of several turbomachines working on their best efficiency point. Plotting these results shows that the best machines lie in a relatively narrow band, so called Cordier or Balje line (cf. figure 2, note that here the definitions of the dimensionless numbers differ slightly from that of equations ( 1) and (2) : σ = 2 -3 /4 π -1 /2 ω s and δ = 2 -3 /4 π 1 /2 D s ). It has to be pointed out that it is a mean empirical line only true in the statistical sense, even if this approach and the resulting classification still a very useful tool for engineers to help guiding the turbomachinery's design from a limited number of inputs. That is to say, it is possible to design a high-efficiency machine far from the Cordier line, especially axial machines as the scatter of points is quite wide in the axial region. Casey et al. [START_REF] Casey | The Cordier Line for Mixed Flow Compressors[END_REF] also outlined several important points in their analysis of the Cordier diagram. As a matter of fact, the dimensionless parameters shown in equations ( 1) and (2) are incompressible by nature. Therefore, the influence of the relative inlet Mach number is not taken into account in the Cordier line. Moreover, the exact location of this line may be strongly driven by the blade geometry, namely the number of blades, the solidity, the pitch angle, the tip clearance and the hub to shroud ratio. Epple et al. [START_REF] Epple | A theoretical derivation of the Cordier diagram for turbomachines[END_REF] derived theoretical expressions for the Cordier line where the influence of the pitch angle and hub to shroud ratio appears clearly and demonstrate the effect of these parameters on the scatter of points on the Cordier diagram. Hence, one could easily consider that the Cordier line and the resulting classification are of less relevance when it comes to design variable geometry or dual operating turbomachines.

A more universal design approach is then needed to handle these new specifications. Lewis [START_REF] Lewis | Turbomachinery Performance Analysis[END_REF] extensively developed a design method based on dimensionless performance variables, namely the work coefficient ψ and the flow coefficient φ defined as:

ψ = ∆h i (ωD) 2 (3)
φ = ṁ ρ i1 ωD 3 (4)
As these parameters and their incompressible derivatives are naturally derived from the dimensional analysis, they often appear in works treating of the unification of turbomachines. For instance, the (ψ, φ) diagram appears in the four-quadrant approaches developed by Knapp et al. [START_REF] Knapp | Complete Characteristics of Centrifugal Pumps and Their Use in the Prediction of Transient Behavior[END_REF] and after that by Turner and Sparkes [START_REF] Turner | Complete characteristics for a single-stage axial flow fan[END_REF] that help exploring all functioning points of a turbomachine (compressor, turbine, dissipation...) whether the flow discharge is positive (design direction) or inversed, and in both rotational directions of the rotors. Gill et al. [START_REF] Gill | Fundamentals of four-quadrant axial flow compressor maps[END_REF] also used the (ψ, φ) diagram through the four-quadrant approach to analyze aerodynamics of compressors under windmilling operation. Binder et al. [START_REF] Binder | Analysis of the Variable Geometry Effect in Radial Turbines[END_REF][START_REF] Binder | Theoretical Analysis of the Aerodynamics of Low-Speed Fans in Free and Load-Controlled Windmilling Operation[END_REF] have been extensively using the (ψ, φ) formalism to predict variable geometry turbine performances and also to treat fan windmilling.

Scope

The (ψ, φ) diagram seems to emerge from the turbomachines' literature when it comes to treat the design of the machine in a unified manner as depicted by Lewis [START_REF] Lewis | Turbomachinery Performance Analysis[END_REF]. In the other hand, this formalism appears as a powerful performance analysis tool regardless of the turbomachine's function or operating point as shown by Gill et al. [START_REF] Gill | Fundamentals of four-quadrant axial flow compressor maps[END_REF] and Binder et al. [START_REF] Binder | Analysis of the Variable Geometry Effect in Radial Turbines[END_REF][START_REF] Binder | Theoretical Analysis of the Aerodynamics of Low-Speed Fans in Free and Load-Controlled Windmilling Operation[END_REF].

Yet, each "type" of turbomachine still, in usual practices of the community, characterized by its own characteristic map that is limited to the narrow windows depicted in figure 1 : The (H, V ) map for hydrodynamic turbomachines, the propeller characteristic map for free or ducted propellers, the (π, ṁ) for compressors and turbines. Therefore, in order to analyze the performances of the current turbomachines in the light of the more unified characterization that represents the (ψ, φ) diagram, it is compulsory to develop relationships between the above mentioned characteristic maps and the (ψ, φ) diagram. The establishment of these relationships represents the main objective of this paper.

Besides, a methodology to extract, from the (ψ, φ) characterization, the turbomachine's flow angles and the meridional plane mean-streamline will be presented. This will allow minimization of the needed inputs for the characteristic maps' conversion relations, hence strengthen their generalization to all turbomachines.

Eventually, the characteristic map conversion and the flow features extraction will be applied to a turbofan's fan stage and also to a propeller which will allow the validation of the developed methods through physical analyses and comparison to experiments. This will also stress the powerful potential of the (ψ, φ) formalism as a performance analysis tool.

Geometric and functioning hypotheses

The focus in this paper will be on single-stage axial compressors and propellers, but the method developed could also be applied to single-stage axial turbines.

Steady, axisymmetric, adiabatic flow of inviscid perfect gas will be considered in the following developments. Furthermore, a mean-line treatment will be made, which means that the layer flow located at a reference radius will be considered as representative of the mean flow in the turbomachine. Therefore, all quantities will be computed for the reference radius defined as [START_REF] Binder | Aéro-thermodynamique des Turbomachines en Fonctionnement Hors-Adaptation[END_REF]:

r j = r 2 sj + r 2 hj 2 (5)
Where subscript s refers to the shroud, h refers to the hub and j refers to the axial location of a given plane (1 for the rotor inlet plane and 2 for the rotor outlet plane).

Regarding the blade-to-blade plane, all angles and velocities of interest are depicted in figure 3. In addition to this, a mean rotor cross-section will be defined as:

S mean = S 1 + S 2 2 (6) W V U W V U U 1 2 1 1 2 2 α1 β2 α2 V θ2 V θ1 1 2 Vz 2 Vz 1 θ z +

Figure 3 Velocity triangles definition

And only on this section, a linear evolution of density and axial velocity will be assumed which leads to:

V zmean = V z1 + V z2 2 (7)
ρ mean = ρ 1 + ρ 2 2 (8)

Characterization maps' definitions

The compressors are often characterized using multiple curves depicting the evolution of the compression ratio and the isentropic efficiency versus corrected mass flow at iso-rotational speed. It will be further referred to this characterization as (π c , M F P, η is ) where M F P is the acronym for Mass Flow Parameter.

π c = p i2 p i1 (compression ratio) (9)
M F P = ṁ S 1 p i1 T i1 R γ ( 10
)
η is = π γ-1 γ c -1 τ c -1 (isentropic efficiency) (11)
It must be stressed here that the rotational speed does not appear in the (π c , M F P, η is ) definition equations. This explains the iso-rotational speed feature of the characteritics in this map and also the necessity of multiple curves to depict the whole operating range of the machine.

The propeller chart characterizes the functioning of a propeller through several curves depicting the evolution of the thrust coefficient and the power coefficient versus the advance ratio at iso-pitch angle. It will be further referred to this characterization as (τ, J, χ).

τ = F ρ i1 U 2 1 S mean (thrust coefficient) (12) J = V z1 U 1 (advance ratio) (13)
χ = P ρ i1 U 3 1 S mean (power coefficient) (14)
Regarding the (ψ, φ) formalism, it depicts the evolution of the work coefficient and the isentropic efficiency versus the flow coefficient at iso-compression ratio. It is then more relevant to refer to this characterization as (ψ, φ, η is ). Moreover the definition further used might slightly differ from that presented in equations ( 3) and ( 4) to be replaced by:

ψ = ∆h i U 2 1 (15) φ = ṁ ρ i1 U 1 S 1 (16)
While the isentropic efficiency η is still defined as in equation [START_REF] Balje | Turbomachines : a guide to design selection and theory[END_REF].

Even if the characterizations previously presented seem different, they share the same physical basis, which is to be derived from the dimensional analysis of a turbomachine. Indeed, the mass flow parameter at iso-rotational speed, the advance ratio or the flow coefficient reflect the rotor inlet relative flow angle, which is with analogy to an airfoil incidence angle, one of the first order driving parameter of the functioning mechanism of the turbomachine. Also, the compression ratio, the thrust coefficient or the work coefficient reflect the quantity of interest that is derived from the machine while the isentropic efficiencies reflect the amount of losses as a turbomachine not only exchanges work with the flow but produces entropy too.

The proximity of these characterizations is exploited in the following developments to relate the different dimensionless parameters by pairs and then create the bridges between the multiple maps.

(π c , M F P, η is ) and (ψ, φ, η is ) maps' relationship Now that all characteristic maps definitions had been set, it appears that passing from the (π c , M F P, η is ) map to the (ψ, φ, η is ) map is kind of straight forward and wellknown. Indeed, having the compressor conventional map and predesign data of geometry and functioning ([i]the rotor inlet radii at hub and shroud; [ii]total pressure, temperature and the flow pre-swirl angle at inlet; [iii]the fluid specific heat ratio) it is possible to express the work coefficient and the flow coefficient as follows:

ψ = c p (T i2 -T i1 ) U 2 1 = γR γ -1 T i1 U 2 1 (τ c -1) (17)
φ = V z1 ρ 1 U 1 ρ i1 (18)
Where the axial velocity V z1 and the static-to-total density ratio ρ1 /ρi 1 at the rotor inlet can be expressed as functions solely depending of the rotor inlet axial Mach number M z1 , through the isentropic flow relations and the Mach number and velocity triangle definitions:

ρ 1 ρ i1 = 1 + γ -1 2 M 2 z1 (1 + tan 2 α 1 ) -1 γ-1 (19)
V z1 = M z1 1 + γ-1 2 M 2 z1 γRT i1 (20)
Eventually the rotor inlet Mach number can be computed using a basic Newton-Raphson method on the following mass conservation equation:

M F P = M z1 1 + γ -1 2 M 2 z1 -γ+1 2(γ-1) (21)
Inversely, by starting from a (ψ, φ, η is ) map and predesign data of geometry and functioning, it is possible to compute the (π c , M F P, η is ) map for a given rotational speed using equations ( 17) to [START_REF] Carter | Low speed performance of related aerofoils in cascades[END_REF].

(τ, J, χ) and (ψ, φ, η is ) maps' relationship

The relationship between the (τ, J, χ) map and the (ψ, φ, η is ) map is more complex than the one previously established. It should then be helpful to first convert the propeller map to a (π c , M F P, η is ) map and then pass it to the (ψ, φ, η is ) map using the previous method.

In order to establish the relation between the (π c , M F P, η is ) map and the (τ, J, χ) map, one must first derive the thrust of the turbomachine from the momentum balance equation:

F = ∂S [(p -p 0 )n + ρ(V • n)V ] dS • z (22)
Then by neglecting the conicality of the system, it is possible to approximate the rotor inlet and outlet sections by the mean rotor cross-section introduced in the equation ( 6), giving:

F ≈ S mean ∆p + ṁ∆V z (23)
Where ∆p = p 2 -p 1 and ∆V z = V z2 -V z1 represents the rotor outlet-inlet pressure and axial velocity variation. This leads, by applying the mass conservation between the rotor inlet, outlet and mean cross-sections and by introducing the equations ( 6), ( 7), ( 8) and ( 13), to:

τ = ∆p ρ i1 V 2 z1 + ρ 1 4ρ i1 (k 2 -1)( 1 ks + 1) J 2 (24)
With k = Vz 2/Vz 1 and s = S2 /S1 are the rotor outlet-to-inlet axial velocity and cross-section ratios respectively. In order to relate the power coefficient χ to the compressor map parameters, one must consider the energy conservation:

P = ṁ∆h i = ṁ γR γ -1 T i1 (τ c -1) (25)
Then by applying the same relations that allow passing from the equation ( 23) to (24) and by using the Mach number definition:

χ = k + 1 γ -1 ρ 1 4ρ i1 1 ks + 1 × (τ c -1) 1 M 2 z1 + γ -1 2 J 3 (26)
Equations ( 13), ( 24) and [START_REF] Mccormick | Aerodynamics, aeronautics, and flight mechanics[END_REF] show that if one knows the (π c , M F P, η is ) map, the rotor outlet section and predesign data of geometry and functioning, it is possible to compute the (τ, J, χ) map. In fact, the advance ratio J can be simply derived by using equations ( 20) and [START_REF] Carter | Low speed performance of related aerofoils in cascades[END_REF]. After that, the determination of the power coefficient χ only needs the computation of the axial velocity ratio k that can be expressed, using the isentropic flow relations and the Mach number definition, as follows:

k = M z2 M z1 τ c 1 + γ-1 2 M 2 z1 1 + γ-1 2 M 2 z2 ( 27
)
Where the rotor outlet axial Mach number is given by applying a Newton-Raphson method on the following mass conservation equation:

M F P √ τ c π c s = M z2 1 + γ -1 2 M 2 z2 -γ+1 2(γ-1) (28)
Still, it is necessary to compute the static pressure variation ∆p to obtain the thrust coefficient τ . This pressure variation can be calculated by using the isentropic flow equations and some basic trigonometric laws applied on the velocity triangles:

p 1 = p i1 1 + γ -1 2 M 2 z1 (1 + tan 2 α 1 ) -γ γ-1 (29)
p 2 = π c p i1 1 + γ -1 2 M 2 z2 × 1 + ( r 2 r 1 1 Jk + tan β 2 ) 2 -γ γ-1 (30)
Inversely, if one starts from a (τ, J, χ) map, the rotor outlet section and predesign data of geometry and functioning, it is possible to compute the (π c , M F P, η is ) map using the previous equations.

Flow angles and the meridional plane meanstreamline extraction

Although some flow features as the relative outlet flow angle β 2 and the mean-streamline rotor outlet radius r 2 are present in the equation (30), they are not mentioned as inputs to pass from a (π c , M F P, η is ) map to a (τ, J, χ) map and vice versa. This is due to the fact that these quantities can be analytically extracted from the (ψ, φ, η is ) map. Indeed, it is possible to relate the work coefficient to the flow coefficient by using the Euler theorem and the velocity triangles definitions as follows:

ψ = ∆h i U 2 1 = U 2 V θ2 -U 1 V θ1 U 2 1 ( 31
)
Which leads, after some manipulations, to:

ψ = r 2 r 1 2 + ρ i1 ρ 1 r 2 r 1 k tan β 2 -tan α 1 φ (32)

ISAIF13

It must be stressed here that the equation ( 32) is not a linear function of the flow coefficient φ since the rotor outletto-inlet axial velocity ratio k will necessarily vary with φ.

The reason is the dependence of k mainly on the compression ratio and the isentropic efficiency as depicted by the equations ( 27) and ( 28). Moreover, one can notice that the factor k will be close to unit in the case of incompressible flow and neglected conicality of the machine, due to mass conservation. The distance of the factor k from unit can therefore be considered as a measure of the incursion of the turbomachine in the compressible domain.

Besides, the rotor outlet relative flow angle β 2 can no more be stated constant if one consider the deviation of the flow from the rotor outlet blade angle. This deviation angle is at iso-geometry of the machine a function of the incidence angle on the blade and Mach number as stated by Carter [START_REF] Carter | Low speed performance of related aerofoils in cascades[END_REF]. This deviation angle can be reformulated as a function of the isentropic efficiency and the k factor, as the former is at first order driven by the incidence angle and the later a function of the variation of the Mach number through the machine.

The non-linearity of the equation (32) can then be related, on the one hand to the compressibility of the flow through the machine and on the other hand to the entropic nature of the flow. The rotor outlet-to-inlet axial velocity ratio k and the isentropic efficiency η is allow to gauge those features of the flow. Consequently, if one consider the optimum operating range of relatively low pressure turbomachines, the variation of k and η is will be small enough to consider the equation (32) as linear.

Moreover, the work coefficient and the flow coefficient can be computed by using exclusively inputs from the rotor inlet plane as shown by the equations ( 17) to [START_REF] Carter | Low speed performance of related aerofoils in cascades[END_REF]. Consequently, the linear fitting of the resulting coefficients and the equation (32) make it possible to extract the values of r 2 and β 2 starting from any characterization map and the limited number of inputs previously enumerated.

Analysis of a converted compressor map

In order to test the coherence of the presented method, it had been applied on the characterization map of the DGEN-380 fan. The DGEN-380 is a turbofan designed to deliver a 2500N thrust and intended to equip personal light jets (4-5 passengers). A complete instrumented test bench of this engine is available at the Department of Aerodynamics, Energetics and Propulsion of ISAE-SUPAERO [START_REF] García Rosa | Experimental Analysis of the Global Performance and the Flow Through a High-Bypass Turbofan in Windmilling Conditions[END_REF]. The figure 4 sums up roughly the architecture of this engine.

The complete fan stage (rotor and stator) had been simulated to produce the conventional (π c , M F P, η is ) map depicted in figure 5. The numerical domain starts upstream of the spinner and ends right upstream of the splitter duct. A mixing plane is located between the rotor and the stator. The geometry of the fan includes the hub fillet and the tip clearance gap.

The numerical domain is discretized with a multi-block structured mesh generated with Numeca's Autogrid 5 T M , w ∼ 1.4 over the blades for all the simulations. The expansion ratio across the boundary layer mesh is about 10%. Non-matching connections are used on periodic faces to facilitate the generation of a high-orthogonality mesh. Overall, the mesh is comprised of about 7.9 million points with 3.5 million points in the rotor and 4.4 million points in the stator.

The simulations are performed with the Euranus solver of the Fine/Turbo T M package of Numeca. For steady problems, the Reynolds-averaged Navier-Stokes equations (RANS) in the rotating frame are solved with a pseudotime-marching method. Time integration is ensured by a four-stage Runge-Kutta scheme with implicit residual smoothing. Local time stepping and a three-level multigrid technique are used to accelerate convergence to the steady state. The discretization in space is based on a cellcentered finite-volume approach. Convective fluxes are determined by a second-order centered scheme with added artificial dissipation of the Jameson type. Viscous fluxes are centered. The mixing plane interface is treated with a conservative coupling by pitchwise rows. Turbulence closure is provided by the one-equation model of Spalart and Allmaras.

At the inlet of the numerical domain, radial profiles are imposed according to the experimental measurements. At outlet, static pressure is imposed with a radial equilibrium condition. For further details about the computational simulation, one can refer to the work of Dufour et al. [START_REF] Dufour | Validation and flow structure analysis in a turbofan stage at windmill[END_REF].

It is interesting to stress here the relatively small compression ratios of this fan in comparison with the conventional aircraft engines' fans. This is due to the smaller diameter of the engine that makes the fan operating in the subsonic domain at all radial stations.

Before applying the map modification, it is crucial to confirm that the starting map is in the scope of the hypotheses made initially. In fact, the presented equations apply only on the rotor of the turbomachine. Before going any further, it must be confirmed that the influence of the stator on the characterization map is negligible. To do so, one must make sure that the losses in the stator are negligible with comparison to the work input of the rotor. Yet these losses are known to rise as the system gets into the choked functioning and make the compression ratio and ISAIF13 the isentropic efficiency plummet. This phenomenon had been observed in the performed numerical simulations and is caused by the formation of shock waves in the stator when the flow is choked. These shock waves raise the losses of the stator in terms of total pressure decrease and make them of the same order of magnitude as the rotor compression ratio. Consequently, by treating exclusively the unchoked part of the compressor map, one can assume the validity of the previously exposed developments as the stator influence on the map can be neglected. First, the unchoked part of the compressor map of the DGEN-380 fan had been converted into a (ψ, φ, η is ) map. The results are provided in the figure 6. The resulting flow and work coefficients are of the same order of magnitude comparing with similar installations [START_REF] Binder | Theoretical Analysis of the Aerodynamics of Low-Speed Fans in Free and Load-Controlled Windmilling Operation[END_REF]. It can be noticed that both ψ and η is characteristics collapses approximatively on the same line as expected. As a matter of fact, this gives the conversion methodology a first level of validation as the uniqueness of these characteristics is predicted by the dimensional analysis theory. As stated before, in the conventional compressor maps, the rotational speed is an independent parameter that is not taken into account in the dimensionless groups as shown by the equations ( 9) to [START_REF] Balje | Turbomachines : a guide to design selection and theory[END_REF], which leads to an isolated characteristic for each rotational speed. Actually, in opposition to the (π c , M F P, η is ) map, the (ψ, φ, η is ) map is the natural result of the dimensional analysis where all independent variables driving the turbomachine's functioning are taken into account except the compressibility and the viscosity. The small discrepancies observed at high flow coefficient are due to these two later parameters that can be represented respectively by the variations of the rotor outlet-toinlet axial velocity ratio k and the isentropic efficiency η is .

Furthermore, the quasi-superposition of the characteristics in the (ψ, φ, η is ) map reflects an operating similarity of the fan at all rotational speeds, which is commonly observed for incompressible hydrodynamic turbomachines. However, the relative tip Mach number of the fan at high rotational speed almost reaches unit value which is, in usual practices of the community, considered as indicating the incursion in the compressible operating domain. Consequently, the relative tip Mach number may not be the relevant parameter to depict the incursion in the compressible domain in the mean-streamline treatment of the turbomachine. The compression ratio and the isentropic efficiency and therefore the factor k seems more relevant to characterize this incursion since they alter the linearity of the equation (32) as stated in the previous section. The table 1 shows that the mean k factor over each iso-rotational speed is close to unit. It is then possible to conclude that the functioning of an unchoked turbomachine can be approximated by an incompressible mean-streamline approach while the variation of the axial velocity through the rotor stays small, even if the relative tip Mach number is high. Given that the (ψ, φ, η is ) map of the DGEN-380 had been computed using the developed method, it make it possible to extract the values of the rotor outlet relative flow angle β 2 and mean-streamline radius r 2 . It is important to stress here that the conversion of the compressor map had been performed at iso-rotational speed. That is to say that the characteristics in the (ψ, φ, η is ) map depict the evolution of the work coefficient and the isentropic efficiency versus the flow coefficient at iso-rotational speed and not at iso-compression ratio. Strictly speaking, the resulting characteristics don't represent a valid (ψ, φ, η is ) map. However, the collapse of the characteristics into a single line due to the little variation of the k factor makes the difference between the iso-rotational speed and the iso-compression ratio characteristics negligible. Still, the isorotational treatment of the map conversion allows building characteristics with a higher number of points which makes the fitting operation far more robust. It is therefore possible to postulate that the higher quality of the fitting operation compensate the physical error introduced by the iso-rotational speed treatment. This is the reason why this little deviation seems acceptable.

The obtained results are very satisfying and give a second level of validation to the method. Actually, the computed flow angle gives β 2 =-32.8 o at the design point of the fan while experiment gives β exp 2 =-29 o [START_REF] García Rosa | Experimental Analysis of the Global Performance and the Flow Through a High-Bypass Turbofan in Windmilling Conditions[END_REF], that is to say a relative error of 13% which is very conclusive for such analytical method. Moreover, the computed meanstreamline radius gives r 2 =136mm that is the exact value obtained by using the equation ( 5), which confirms the relevance of this definition of the mean-streamline.

The computation of the values of the rotor outlet relative flow angle and mean-streamline radius makes it possible to convert the initial (π c , M F P, η is ) map to (τ, J, χ) map by applying the presented method. The obtained results are depicted in the figure 7. One can notice that the resulting thrust coefficient, power coefficient and advance ratio are of the same order of magnitude comparing with propellers producing equivalent thrust (i.e. the 568F propeller driven by the PW-127M engine that equips the ATR-72 [24]). Similarly to the observation made on the (ψ, φ, η is ) map, the characteristics on the propeller map collapse into a single line. This result confirms the validity of the developed method since the (τ, J, χ) map takes also into account all the independent variables driving the turbomachine's functioning except the compressibility and the viscosity. Also, this quasi-superposition of the characteristics confirms the small influence of the compressibility on the global functioning parameters of the fan as analyzed before.

Analysis of a converted propeller map

In the previous section, the developed methodology had been tested by converting a compressor map into (ψ, φ, η is ) map and propeller map. Still, it is interesting to test the bijectivity of the method by starting from a propeller map. To do so, the wind tunnel test results of the 3.048m diameter four blade Clark-Y 5868-9 propeller performed by Hartman and Biermann [START_REF] Hartman | The aerodynamic characteristics of full-scale propellers having 2, 3, and 4 blades of clark Y and R.A.F. 6 airfoil sections[END_REF] had been taken as reference propeller map. Hartman and Biermann produced the complete propeller characterization for multiple blade pitch angles and presented it in a classical fashion through a (τ, J, χ) map that is depicted in figure 8. The conversion methodology had been applied to each iso-pitch angle characteristic and the results are depicted in the figure 9. Here, the mean k factor value on each curve is about 0.98 which reflects the incompressible aspect of the flow through the propeller. That allows admitting the validity of the map conversion along iso-pitch angle curves. First it can be noticed that the propeller iso-pitch angle characteristics in the (ψ, φ, η is ) map are linear as expected by the theoretical derivation shown in equation (32). In this case, the uniqueness of the characteristics is not observed due to the variation of the pitch angle that represents an independent variable not taken into account in the dimensional analysis performed to build the dimensionless ISAIF13 groups. Also, even if the resulting characteristics seem to be parallel, it is actually not the case. The table 2 sums up the linear fitting coefficients of each characteristic and shows that the characteristics are not parallel. It is interesting to focus on the coefficient b that depicts, according to the equation (32), the squared rotor outlet-toinlet mean-streamline radius ratio. It takes under unit values which can seem unexpected. Indeed, in compressors for which the (ψ, φ, η is ) map is more often used, the convergence of the cross-section is a design commonly admitted as it allows stable operating of the machine by maintaining a quasi-constant axial velocity. This feature, in addition to the casing that confines the flow, makes the meanstreamline increase radius which leads to a rotor outletto-inlet mean-streamline radius ratio over unit. However, in the case of the studied propeller, there is neither crosssection evolution nor casing which makes the stream tube's geometry defining the flow through the propeller function of its operating parameters. As a matter of fact, this stream tube has a convergent form due to the propeller work input [START_REF] Mccormick | Aerodynamics, aeronautics, and flight mechanics[END_REF]. This is the reason why the computed rotor outlet-toinlet mean-streamline radius ratio takes values under unit. Moreover, as depicted in the reference literature [START_REF] Mccormick | Aerodynamics, aeronautics, and flight mechanics[END_REF], the convergence of the stream tube decreases with the increase of the advance ratio. This phenomenon can also be observed in the table 2 through the evolution of the value of the rotor outlet-to-inlet mean-streamline radius ratio which tends to unit when the pitch angle and therefore the related advance ratio increase (cf. figure 8).

These observations demonstrate that the meanstreamline radius definition given by the equation ( 5) is not relevant for the treatment of free propellers; one better prefers in that case simply the rotor radius as it is representative of the stream tube radius at upstream far field.

Eventually, the rotor outlet relative flow angle β 2 had been extracted using the developed method and had been compared to the blade pitch angles considering that this angle reflects the mean blade profile trailing edge angle. The table 3 sums up the obtained results. It is important to clarify that these results are not meant to be quantitative as more experimental data are needed to make a proper comparison. The objective here is to stress the ability of the (ψ, φ, η is ) formalism to highlight off-design functioning.

As expected, except for the 30 o pitch angle which is the design pitch of the propeller, the values of the deviation angle are quite high. Indeed, the low solidity of the propellers due to their little number of blades make these machines far more sensitive to flow detachments when operating at off-design points.

The consistency of the presented observations with the well-known physical phenomena allows stressing the validity of the developed map conversion methods. Also, one must focus on the (ψ, φ, η is ) map as a powerful turbomachines' performance analysis tool. Actually, it is that formalism that makes it possible to quantify the compressibility effect, to extract kinematic information about the flow field -therefore about the turbomachine's geometry -and to enlighten off-design operating points.

Conclusions

The present study had undertaken the definition of analytical relationships between the different turbomachines' characterization maps, namely the compressor/turbine (π c , M F P, η is ) map, the propeller (τ, J, χ) map and the (ψ, φ, η is ) map. The developed equations are bijective which makes it possible to convert any initial characteristic map into the other two in order to take advantage of each map. Actually, the (ψ, φ, η is ) map had been used to develop a mean flow features extraction method that determines the rotor outlet relative flow angle and the meanstreamline radius.

These relationships and the mean flow features extraction method had been applied to both turbofan's fan stage and a propeller which helped to highlight the powerful potential of the (ψ, φ, η is ) formalism as a performance analysis tool. Indeed, the physical analyses of the obtained results and comparisons to experiments allow drawing some important conclusions:

• The compressibility effects on the global performance parameters of an unchoked turbomachine are negligible while the variation of axial flow stays small (i.e. k factor close to unit).

• The rotor outlet relative flow angle at reference radius of a single-stage unchoked fan can be analytically derived from its characteristic map with an error up to 15%.

• The definition of the meridional mean-streamline of the flow in a single-stage unchoked fan can be obtained through the equation ( 5) as proposed by Binder [START_REF] Binder | Aéro-thermodynamique des Turbomachines en Fonctionnement Hors-Adaptation[END_REF].

• The geometry of the stream tube that defines the flow through a propeller can be analytically derived using the equation (32), both for design and off-design operating points.

• The design pitch angle of a propeller can be analytically derived through the extraction of the flow angles and the analysis of the resulting deviation angle.

These results, through their coherence with the physics of the studied turbomachines bring a first level of validation to the developed equations. Finally, the developed methods allow the characterization of all axial turbomachines using the most unified tool that is the (ψ, φ, η is ) formalism. This map already showed its relevance in the unification of the turbomachines' treatment through its usage in the works of Gill et al. [START_REF] Gill | Fundamentals of four-quadrant axial flow compressor maps[END_REF] and Binder et al. [START_REF] Binder | Analysis of the Variable Geometry Effect in Radial Turbines[END_REF][START_REF] Binder | Theoretical Analysis of the Aerodynamics of Low-Speed Fans in Free and Load-Controlled Windmilling Operation[END_REF][START_REF] Binder | Aéro-thermodynamique des Turbomachines en Fonctionnement Hors-Adaptation[END_REF] on low-speed fans and radial turbines. This contribution adds turbofan's fan stages and propellers to the later list. In addition to this, it is important to stress the ability of the (ψ, φ, η is ) map to represent clearly the off-design operation of the turbomachine. In particular the blades pitch angle variation that opens very promising prospects for the analytical modelisation of this geometrical variability inside current system-oriented reasearch on innovative propulsive systems.

  Nomenclature c p : specific heat at constant pressure [J/kg/K] D : rotor diameter [m] D s : specific diameter [-] F : thrust [N] J : advance ratio [-] k : rotor axial velocity ratio [-] ṁ : mass flow rate [kg/s] n : conventional unit normal vector [-] N : rotational speed [RPM] p : pressure [Pa] P : power [W] r : rotor radial station [m] R : gas constant [J/kg/K] s : rotor outlet-to-inlet cross-section ratio [-] S : cross-section [m 2 ] T : temperature [K] U : tangential rotational velocity [m/s] V : absolute velocity [m/s] V : flow rate [m 3 /s] z : unit axial vector [-] α : absolute flow angle [rad] β : relative flow angle [rad] γ : specific heat ratio [-] ∆h i : total enthalpy variation [J] η : efficiency [-] π c : total-to-total compression ratio [-] ρ : density [kg/m 3 ] τ : thrust coefficient [-] τ c : total-to-total temperature ratio [-] φ : flow coefficient [-] χ : power coefficient [-] ψ : work coefficient [-] ω : rotational speed [rad/s] ω s : specific speed [-] Subscripts h : hub i : total quantity is : isentropic j : axial station nom : nominal opt : optimal s : shroud z : axial projection

Figure 1

 1 Figure 1 Turbomachines design windows (adapted from Japikse & Baines [6])

Figure 2

 2 Figure 2 Cordier diagram [10]

Figure 4

 4 Figure 4 DGEN-380 architecture consisting of O-type blocks around the blades, and H-type blocks to fill the passage. The first cell at the wall has a size of 10 µm, ensuring a value of ∆y +w ∼ 1.4 over the blades for all the simulations. The expansion ratio across the boundary layer mesh is about 10%. Non-matching connections are used on periodic faces to facilitate the generation of a high-orthogonality mesh. Overall, the mesh is comprised of about 7.9 million points with 3.5 million points in the rotor and 4.4 million points in the stator.The simulations are performed with the Euranus solver of the Fine/Turbo T M package of Numeca. For steady problems, the Reynolds-averaged Navier-Stokes equations (RANS) in the rotating frame are solved with a pseudotime-marching method. Time integration is ensured by a four-stage Runge-Kutta scheme with implicit residual smoothing. Local time stepping and a three-level multigrid technique are used to accelerate convergence to the steady state. The discretization in space is based on a cellcentered finite-volume approach. Convective fluxes are determined by a second-order centered scheme with added artificial dissipation of the Jameson type. Viscous fluxes are centered. The mixing plane interface is treated with a conservative coupling by pitchwise rows. Turbulence closure is provided by the one-equation model of Spalart and Allmaras.At the inlet of the numerical domain, radial profiles are imposed according to the experimental measurements. At outlet, static pressure is imposed with a radial equilibrium condition. For further details about the computational simulation, one can refer to the work of Dufour et al.[START_REF] Dufour | Validation and flow structure analysis in a turbofan stage at windmill[END_REF].It is interesting to stress here the relatively small compression ratios of this fan in comparison with the conventional aircraft engines' fans. This is due to the smaller diameter of the engine that makes the fan operating in the subsonic domain at all radial stations.Before applying the map modification, it is crucial to confirm that the starting map is in the scope of the hypotheses made initially. In fact, the presented equations apply only on the rotor of the turbomachine. Before going any further, it must be confirmed that the influence of the stator on the characterization map is negligible. To do so, one must make sure that the losses in the stator are negligible with comparison to the work input of the rotor. Yet these losses are known to rise as the system gets into the choked functioning and make the compression ratio and

Figure 5

 5 Figure 5 DGEN-380 (πc, M F P, η is ) map

Figure 6

 6 Figure 6 DGEN-380 (ψ, φ, η is ) map

Figure 7

 7 Figure 7 DGEN-380 (τ, J, χ) map

Figure 8

 8 Figure 8 Clark-Y 5868-9 (τ, J, χ) map

Figure 9

 9 Figure 9 Clark-Y 5868-9 (ψ, φ, η is ) map

Table 1

 1

	DGEN-380 rotor mean outlet-to-inlet axial velocity
	ratio k
	Rotational speed [%Nnom]	k
	60	1.0712
	70	1.0631
	90	1.0165
	95	1.0012
	100	0.9860
	108	0.9612
	118	0.9245

Table 2

 2 Clark-Y 5868-9 (ψ, φ, η is ) map linear fitting coefficients as ψ = aφ + b

	Pitch angle	a	b
	20 o	-0.4635 0.4404
	25 o	-0.4657 0.5272
	30 o	-0.4740 0.6386
	35 o	-0.4566 0.7192
	40 o	-0.5601 0.9741

Table 3

 3 Flow and blade angle comparison

	Pitch angle Flow angle Deviation angle
	20 o	33.8 o	13.8 o
	25 o	30.9 o	5.9 o
	30 o	28.2 o	1.8 o
	35 o	25.6 o	9.4 o
	40 o	25.4 o	14.6 o


Alexander W Winkler
email: alexander.winkler@mavt.ethz.ch

Carlos Mastalli
email: carlos.mastalli@iit.it.

Ioannis Havoutis
email: ioannis.havoutis@iit.it.

Michele Focchi
email: michele.focchi@iit.it.

Darwin G Caldwell
email: darwin.caldwell@iit.it.

Claudio Semini
email: claudio.semini@iit.it.

Planning and Execution of Dynamic Whole-Body Locomotion for a Hydraulic Quadruped on Challenging Terrain

We present a framework for dynamic quadrupedal locomotion over challenging terrain, where the choice of appropriate footholds is crucial for the success of the behaviour. We build a model of the environment on-line and on-board using an efficient occupancy grid representation. We use Any-time-Repairing A* (ARA*) to search over a tree of possible actions, choose a rough body path and select the locally-best footholds accordingly. We run a n-step lookahead optimization of the body trajectory using a dynamic stability metric, the Zero Moment Point (ZMP), that generates natural dynamic whole-body motions. A combination of floating-base inverse dynamics and virtual model control accurately executes the desired motions on an actively compliant system. Experimental trials show that this framework allows us to traverse terrains at nearly 6 times the speed of our previous work, evaluated over the same set of trials.

I. INTRODUCTION

Agile locomotion is one of the key abilities that legged ground robots need to master. Wheeled or tracked vehicles are efficient in structured environments but can suffer from limited mobility in many real-world scenarios. Legged robots offer a clear advantage in unstructured and challenging terrain. Such environments are common in disaster relief, search & rescue, forestry and construction site scenarios.

This paper presents the newest development in a stream of research that aims to increase the autonomy and flexibility of legged robots in unstructured and challenging environments. We present a framework for dynamic quadrupedal locomotion over highly challenging terrain where the choice of appropriate footholds is crucial for the success of the behaviour. We use perception to build a map of the environment, decide on a rough body path and choose appropriate footholds. We are able to generate feasible footholds on-line and on-board for various types of scenarios such as climbing up and down pallets, traversing stepping stones using an irregular swingleg sequence and passing over a 35 cm gap. We optimize the body trajectory according to a dynamic stability metric (ZMP) to produce agile and natural dynamic whole-body motions up to 5.8 times the speed of our previous work [START_REF] Winkler | Path planning with force-based foothold adaptation and virtual model control for torque controlled quadruped robots[END_REF]. Compliant execution of the motions is performed using a floating-base inverse dynamics controller that ensures the accurate execution of dynamic motions, in combination with a virtual model controller that generates feedback torques to account for model and tracking inaccuracies. The hydraulically actuated and fully torque controlled quadruped robot HyQ [START_REF] Semini | Design of HyQ -a hydraulically and electrically actuated quadruped robot[END_REF]. The inset plot shows the on-line built environment perception alongside with the planned footholds and the current on-board robot state estimate.

Our contribution includes on-line perception, map building and foothold planning, generation and execution of optimized dynamic whole-body motions despite irregular swing-leg sequences and the use of an elegant inverse-dynamics/virtualmodel control formulation that exploits the natural partitioning of the robot's dynamic equations.

The rest of the paper is structured as follows: After discussing previous research in the field of dynamic wholebody locomotion (II) we describe the on-line map building and how appropriate footholds are chosen (III). Section IV explains, how dynamically stable whole-body motions are generated based on an arbitrary footstep sequence. Section V shows how these desired motions are accurately and compliantly executed. In Section VI we evaluate the performance of our framework on the Hydraulic Quadruped robot HyQ (Fig. 1) in real-world experimental trials before Section VII summarizes this work and presents ideas for future work.

II. RELATED WORK

In environments where smooth, continuous support is available (flats, fields, roads, etc.), where exact foot placement is not crucial for the success of the behaviour, legged systems can utilize a variety of dynamic gaits, e.g. trotting, galloping. Marc Raibert pioneered the study of the principles of dynamic balancing with legged robots [START_REF] Raibert | Legged robots that balance[END_REF], resulting in the quadruped BigDog. The reactive controllers used in these legged systems are partially capable to overcome unstructured terrain. Likewise, HyQ can traverse lightly unstructured terrain using reactive control [START_REF] Barasuol | A reactive controller framework for quadrupedal locomotion on challenging terrain[END_REF], [START_REF] Havoutis | Quadrupedal trotting with active compliance[END_REF] or reflex strategies [START_REF] Focchi | Local Reflex Generation for Obstacle Negotiation in Quadrupedal Locomotion[END_REF].

However, for more complex environments with obstacles like large gaps or stairs, such systems quickly reach their limits. In this case, higher level motion planning that considers the environment and carefully selects appropriate footholds is required. In these terrains, e.g. stairs, gaps, cluttered rooms, legged robots have the potential to use non-gaited locomotion strategies that rely more on accurate foothold planning based on features of the terrain. There exist a number of successful control architectures [START_REF] Kalakrishnan | Learning, planning, and control for quadruped locomotion over challenging terrain[END_REF], [START_REF] Rebula | A controller for the littledog quadruped walking on rough terrain[END_REF], [START_REF] Kolter | A control architecture for quadruped locomotion over rough terrain[END_REF], [START_REF] Zucker | Optimization and learning for rough terrain legged locomotion[END_REF] to plan and execute footsteps to traverse such terrain. Some avoid global footstep planning by simply choosing the next best reachable foothold [START_REF] Rebula | A controller for the littledog quadruped walking on rough terrain[END_REF], while others plan the complete footstep sequence from start to goal [START_REF] Zucker | Optimization and learning for rough terrain legged locomotion[END_REF], often requiring time consuming re-planning in case of slippage or deviation from the planned path.

The approach in [START_REF] Kalakrishnan | Learning, planning, and control for quadruped locomotion over challenging terrain[END_REF] stands between the two above mentioned methods and plans a global rough body path to avoid local minima, but the specific footholds are chosen only a few steps in advance. This reduces the necessary time for replanning in case of slippage, while still considering a locally optimal plan. We recently built on this approach with a path planning and control framework that uses on-line force-based foothold adaptation to update the planned motion according to the perceived state of the environment during execution [START_REF] Winkler | Path planning with force-based foothold adaptation and virtual model control for torque controlled quadruped robots[END_REF].

The whole-body locomotion framework described in this paper further extends this work: We use real-time perception to create, evaluate and update a terrain cost map on-board. Compared to previous approaches our framework does not make use of any external state measuring system, e.g. a marker-based tracking system. The incorporation of domain knowledge, e.g. body motion primitives and an ARA* planner, allows us to re-plan actions and footholds on-line. As in [START_REF] Kalakrishnan | Learning, planning, and control for quadruped locomotion over challenging terrain[END_REF] the Center of Gravity (CoG) trajectory is now chosen to comply with the ZMP dynamic stability metric [START_REF] Kajita | Biped walking pattern generation by using preview control of zero-moment point[END_REF] to produce agile, fast and natural motions.

III. PERCEPTION AND (RE-)PLANNING

This section describes the pipeline from the acquisition and evaluation of terrain information to the generation of appropriate footholds (Fig. 2) 1 . The on-board terrain information server continuously holds the state of the environment. The body action planner decides the general direction of movement and the footstep sequence planner chooses specific footholds along this path.

A. Terrain Information

We develop a terrain information server that computes the required information for the body action and the footstep sequence planners, e.g. the terrain cost map of the environment. We build a 3D occupancy grid map [START_REF] Hornung | OctoMap: An efficient probabilistic 3D mapping framework based on octrees[END_REF] from a RGBD sensor mounted on a scanning pan & tilt unit, alongside with the state estimate of the robot using the Extended Kalman Filter [START_REF] Bloesch | State estimation for legged robotsconsistent fusion of leg kinematics and imu[END_REF]. The voxel-based map is built using a (x, y, z) resolution of (4 cm × 4 cm × 2 cm) which roughly matches the dimensions of the robot's foot. 1 A more in-depth presentation of the perception and terrain evaluation pipeline can be found in [START_REF] Mastalli | On-line and on-board planning for quadrupedal locomotion, using practical, on-board perception[END_REF]. The terrain cost map quantifies how desirable it is to place a foot at a specific location. The terrain cost c t for each voxel in the map is computed using geometric terrain features as in [START_REF] Winkler | Path planning with force-based foothold adaptation and virtual model control for torque controlled quadruped robots[END_REF]. Namely, we use the standard deviation of height values, the slope and the curvature of the cell in question. The terrain cost c t for each voxel of the map is computed as a weighted linear combination of the individual features c t (x, y) = w T c(x, y). The cost map is locally re-computed (in a 2.5 m×5 m area around the robot) whenever a change in the map is detected.

B. Body Action Planner

The state s = (x, y, θ) ∈ S of the robot body includes the current position (x, y) and the yaw angle θ. Given a desired goal state, the body action planner finds a sequence of actions a 0...N = {a 0 , a 1 , . . . , a N } that move the robot in a nearly optimal way to this state. This implies that terrain features, the difficulty of specific actions, kinematic reachability and collision with the environment must be considered and quantified. A feasible action a is the change of state that can be achieved through one step from s to s as a = (∆x, ∆y, ∆θ) ∈ A.

We define A as the set of -empirically chosen-feasible motion primitives (e.g. move left, diagonally forward, back) that correspond to the kinematics and dynamics of the robot. The cost of an action a given a current state s is computed as a weighted linear combination of costs:

c(s, a) = w T c(s, a) (2)
with c(s, a) consisting of: ct The average of the best n terrain costs around each leg after performing action a. c a

The difficulty of a specific action, e.g. sideways steps are more difficult than forward ones. c pc Penalizes actions that potentially cause the swingleg to collide with the environment. c po Penalizes actions that potentially end up in uneven terrain that require large roll and pitch angles. The objective is to find a sequence of actions a from the current body state s = (x, y, θ) to the goal state g, that minimizes the accumulated action costs c(s, a). For simplicity only three possible actions are shown, namely move left (a l ), right (a r ) and forward (a f ). The optimal action sequence {a l , a r , . . . , r } found through ARA* is shown in red.

The set of actions A and the current state s of the robot is used to construct a directed graph G = (S, A) (Fig. 3). We use the ARA* [START_REF] Likhachev | Ara*: Anytime a* with provable bounds on sub-optimality[END_REF] algorithm to search the tree for a sequence of actions with the lowest accumulated cost from the current to the goal state. ARA* uses a heuristic h(s) = -cF( gs ) to decide along which states to search first. F(•) is the estimated remaining steps (actions) to reach the goal state and c is an estimated lower bound on the average future action costs, considering the terrain costs between the current and the goal state.

ARA* initially runs an A* search with an inflated heuristic, • h(s), which quickly finds a first sequence of actions. Unfortunately, since the inflated heuristic is no longer admissible (always lower than the true cost), the sequence of actions may be sub-optimal. As long as computational time is still available, ARA* repeatedly runs A* search, continuously decreasing the inflation factor and thereby finding closer to optimal sequences of actions. Since a first solution, although suboptimal, is found quickly, this algorithm can be used online.

C. Footstep Sequence Planner

Given the desired body action plan, the footstep sequence planner computes the sequence of footholds that corresponds to these body actions. In our previous work, we selected the optimal foothold around the nominal stance positions in a predefined swing-leg sequence. In this paper, we modify the position of the search area and the swing-leg sequence depending on the corresponding action, which improves the robustness of the planned actions. For example, when moving left (action a l ) it is advantageous to swing one of the left legs to avoid small areas of support.

The footstep location in each search area with the lowest foothold cost, c f = w t c t + w st c st + w c c c + w o c o , is then selected, where c t is the terrain cost below a foothold, c st is the support triangle cost, c c is the leg collision cost and c o is the body orientation cost.

D. Re-planning/Updating during Execution

Compared to our previous approach the graph is significantly smaller, since we only search over feasible actions A and not over every discretized change in state. Additionally, ARA* provides intermediate solutions, so the exhaustive and time-costly search procedure does not need to be completed before the robot can react. This combination of the efficient voxel-based occupancy map, the graph representation over feasible actions, and the efficient search through ARA* allows us to re-plan the motions and set of planned footholds online to cope with changes in the environment as is shown in Fig. 4.

IV. WHOLE-BODY MOTION GENERATION

We generate a body trajectory that ensures that the robot is dynamically stable at every time step. We follow the approach presented in [START_REF] Kalakrishnan | Learning, planning, and control for quadruped locomotion over challenging terrain[END_REF] that finds a CoG trajectory that respects stability constraints without explicitly generating a ZMP trajectory. We build on this approach by enabling swing-leg sequences in any order through insertion of fourleg support phases.

A. Problem Formulation

For a CoG trajectory to be feasible it must be continuous and double differentiable. This way we avoid steps in accelerations that produce discontinuous torques which can damage the hardware and affect stability. The body trajectory, x cog , is given by a spline comprised of multiple fifth-order polynomials:

x cog (t) = a x t 5 + b x t 4 + c x t 3 + d x t 2 + e x t + f x . (3)
At each spline junction we require the last state (t = T i ) of spline i to be equal to the first state (t = 0) of the next spline i+1 as:

(x cog , ẋcog , ẍcog ) i t=Ti = (x cog , ẋcog , ẍcog ) i+1 t=0 . (4)
This ensures double differentiability and continuity of the trajectory, required by the floating-base inverse dynamics. Finding an optimal CoG trajectory can then be reduced to finding optimal polynomial coefficients q i = (a x , . . . , f x , a y , . . . , f y ) T ∈ R 12 for each spline segment i.

B. Dynamic Stability

To execute the planned footsteps, a body trajectory must be found that ensures a stable stance at all time instances. For slow movements this is achieved by keeping the CoG inside the support polygon, i.e. the polygon formed by the legs in stance. To consider dynamic effects of larger body accelerations we estimate the position of the ZMP by modeling the robot as a cart-table (Fig. 5, left). The ZMP can then be calculated by:

x zmp = x cog - z cog ẍcog zcog + g 0 , (5)
where x zmp and x cog are the position of the ZMP and the CoG respectively, z cog describes the height of the robot with respect to its feet, zcog is the vertical acceleration of the body and g 0 represents the gravitational acceleration. Dynamic stability requires the ZMP to be inside the current support triangle, expressed by three lines l of the form px + qy + r = 0. The ZMP is considered to be inside a support triangle, if the following conditions are met at every sampling interval: p l x zmp + q l y zmp + r l > 0 for l = 1, 2, 3.

In reality there exist discrepancies between the cart-table model and the real robot. Additionally, desired body trajectories cannot be perfectly tracked as modelling, sensing and actuation inaccuracies are hard to avoid. Therefore, it is best to avoid the border of stable configurations by shrinking the support triangles by a stability margin d (Fig. 5, right). With the introduction of d, there is no continuous ZMP trajectory when switching between diagonally opposite swing legs as the support triangles are disjoint. We therefore allow a transition period ('four-leg support phase') during which the ZMP is only restricted by the shrunk support polygon created by the four stance feet.

We built on [START_REF] Kalakrishnan | Learning, planning, and control for quadruped locomotion over challenging terrain[END_REF] by allowing a completely irregular sequence of steps for the ZMP optimization. Our trajectory generator needs no knowledge of a predefined gait. For every step it checks if the next swing leg is diagonally opposite of the current swing leg. If so, the disjoint support triangles require a four-leg support phase for the optimization to find a solution. This allows a greater decoupling from the footstep planner, which can generate swing leg sequences in any order useful for the success of the behaviour.

C. Cost function

In addition to moving in a dynamically stable way, the trajectory should accelerate as little as possible during the execution period T . This increases possible execution speed and reduces required joint torques. This is achieved by minimizing

J = w x T 0 ẍ2 cog (t) dt + w y T 0 ÿ2 cog (t) dt. (7)
The directional weights w penalize sideways accelerations (w y = 1.5w x ) more than forward-backward motions, since sideways motions are more likely to cause instability. This results in a convex quadratic program (QP) with the cost function [START_REF] Kalakrishnan | Learning, planning, and control for quadruped locomotion over challenging terrain[END_REF], the equality constraints (4) and the inequality constraints [START_REF] Focchi | Local Reflex Generation for Obstacle Negotiation in Quadrupedal Locomotion[END_REF]. We solve it using the freely available QP solver, namely QuadProg++ [START_REF] Guennebaud | eiquadprog.hh[END_REF] to obtain the spline coefficients q and therefore the desired and stable (x, y)-body trajectory (3). The remaining degrees of freedom (z cog , roll, pitch, yaw) and the swing-leg trajectories are chosen based mostly on the foothold heights and described in detail in [START_REF] Winkler | Path planning with force-based foothold adaptation and virtual model control for torque controlled quadruped robots[END_REF].

V. EXECUTION OF WHOLE-BODY MOTIONS

Dynamic whole-body motions require orchestrated and precise actuation of all the joints. Simple PD controllers do not suffice for such motions, especially when considering uncertainties in the environment and/or model inaccuracies. We use a control scheme (Fig. 6) that combines a virtual model with a floating-base inverse dynamics controller. After receiving an arbitrary sequence of footholds from the footstep planner, the whole-body motion generator calculates desired (feedforward) accelerations ẍd for the body and a virtual model (VM) control loop adds feedback accelerations ẍfb should the robot deviate from the desired trajectory. The inverse dynamics produce the majority of the joint torques which are combined with a low-gain joint-space PD controller to compensate for possible model inaccuracies. The computed reference torques are then tracked by the lowlevel torque controller. Note that x describes the linear and rotational coordinates of the body as [START_REF] Rebula | A controller for the littledog quadruped walking on rough terrain[END_REF] where R b ∈ R 3×3 is a coordinate rotation matrix representing the orientation of the base w.r.t. the world frame and ω b ∈ R 3 is the angular velocity of the base.

x = (x cog , R b ), ẋ = ( ẋcog , ω b ), ẍ = (ẍ cog , ωb ),

A. Virtual Model

The feedback action which compensates for inaccurate execution and drift can be imagined as virtual springs and dampers attached to the robot's trunk on one side and the desired body trajectory on the other [START_REF] Pratt | Virtual model control: An intuitive approach for bipedal locomotion[END_REF]. Deviation between these causes the springs and dampers to produce virtual forces F vm and torques T vm on the body that "pull" the robot back into the desired state through

F vm = P x (x d cog -x cog ) + D x ( ẋd cog -ẋcog ) T vm = P θ e(R d b R b ) + D θ (ω d b -ω b ), (9)
where the superscript d refers to the desired values, planned by the whole-body motion generator and non-superscript values describe the state of the robot as estimated by the on-board state estimator. Respectively, e(.) : R 3×3 → R 3 is a mapping from a rotation matrix to the associated rotation vector [START_REF] Caccavale | Six-DOF impedance control based on angle/axis representations[END_REF]. P x , D x , P θ , D θ ∈ R 3×3 are positive-definite diagonal matrices of proportional and derivative gains, respectively. Expressing the body feedback action in terms of forces and moments allows us give the virtual model gains a physical meaning of stiffness and damping and thus can be intuitively set and used. Since the inverse dynamics computation requires reference accelerations, we multiply forces/moments (wrench) W vm = (F vm , T vm ) with the inverse of the composite rigid body inertia I c of the robot at its current configuration. Adding this body feedback acceleration to the desired body acceleration produced by the whole-body motion generator creates the 6D reference acceleration (linear and rotational) for the inverse dynamics computation as:

ẍref = ẍd + I -1 c W vm . (10)
By combining a feedforward acceleration ẍd with a bodyfeedback acceleration, we achieve accurate tracking while maintaining a compliant behaviour.

B. Floating Base Inverse Dynamics

The floating base inverse dynamics algorithm calculates the joint torques required to execute the reference body accelerations. We can partition [START_REF] Fujimoto | Robust biped walking with active interaction control between foot and ground[END_REF] the dynamics equation of the robot into the unactuated base coordinates q b ∈ R 6 and the active joints' q ∈ R 12 as

M(R, q) qb q + h b h q (R, q, ω, q) b = 0 τ + J T cb J T cq λ, (11)
where M is the floating base mass matrix, h = (h b , h q ) is the force vector that accounts for Coriolis, centrifugal, and gravitational forces, λ are the ground contact forces, and their corresponding Jacobian J c = J cb J cq and τ are the torques that we wish to calculate.

The left hand term b = (b b , b q ) can be computed efficiently using the Featherstone implementation of the Recursive Newton-Euler Algorithm (RNEA) [START_REF] Featherstone | Rigid Body Dynamics Algorithms[END_REF]. Since the CoG acceleration ẍref cog is defined in a frame aligned with the base frame but with the origin in the CoG, we perform a translational coordinate transform b X cog to get the 6D base spatial acceleration: qb = b X cog ẍref as in [START_REF] Featherstone | Rigid Body Dynamics Algorithms[END_REF]. By partitioning the dynamics equation as in [START_REF] Kajita | Biped walking pattern generation by using preview control of zero-moment point[END_REF] and given that the base is not actuated, we can directly compute, in a least-squares way, the vector of ground reaction forces λ from the first n b equations, λ = J + cb b b , where () + denotes the Moore-Penrose generalized inverse. We then use the last n equations to produce the reference joint torques, τ id = b q -J T cq λ.

VI. EXPERIMENTAL RESULTS

This section describes the experiments conducted to validate the effectiveness and quantify the performance of our framework.

A. Experimental Setup

We use the hydraulically-actuated quadruped robot HyQ in our experiments. HyQ weighs approximately 90 kg, is fullytorque controlled and equipped with precision joint encoders, a depth camera (Asus Xtion) and an Inertial Measurement Unit (MicroStrain). We perform on-board state estimation and do not make use of any external state measuring system, e.g. a marker-based tracking system. All computations are done on-board, using a PC104 stack for the real-time critical part of the framework, and a commodity i7/2.8 GHz PC for perception and planning.

The first experiment starts with flat, obstacle-free terrain. After the robot has planned initial footsteps, a pallet is placed into the terrain. In the next experiments the robot must climb one and two pallets of dimensions 1.2 m × 0.8 m × 0.15 m. The height of one pallet is equal to 20% of the leg length. Furthermore we show that the robot traverses a gap of 35 cm, which is approximately half the robot's body length. The final experiment consist of two pallets connected by a sparse path of stepping stones. The pallets are 1.2 m apart and the stepping stones lie 0.08 m lower than the pallets.

For each experiment, we specify the (x, y, θ) goal state. The footstep planner finds a sequence of footsteps of an arbitrary order, which the controller then executes dynamically. We validate the performance of our framework in 4 scenarios as seen in Fig. 7 and compare it to our previously achieved results (Table I) on the same benchmark tasks. Additionally, the reader is strongly encouraged to view the accompanying video 2 as it provides the most intuitive way to demonstrate the performance of our framework.

B. Results and Discussion

1) Perception and (re-)planning: Efficient occupancy grid-based mapping and focusing our computations to an area of interest around the robot body greatly increase computation speed. This allows us to incrementally build a model of the environment and update the terrain cost map at a frequency of 2 Hz. Using the action based search I) that our framework can achieve: This is due to the fact, that the body can move faster while still being stable, since we are using a dynamic stability criterion. All accelerations and decelerations are optimized, so that the ZMP never leaves the support polygon. In addition, since we are not directly producing torques with the virtual model feedback controller, but only accelerations for the inverse dynamics controller, our feedback actions also respect the dynamics of the system. Furthermore, the duration of the four-leg-support phase is significantly reduced: It is much faster to move the ZMP from one support triangle to another than the CoG (e.g. entire body), because this can be achieved by manipulating the acceleration.

3) Model accuracy: Walking over a 35 cm gap (approximately half of the body length) shows the stability of the robot despite of highly dynamic motions. When crossing the gap the robot accelerates up to a body velocity of 0.5 m/s and is able to decelerate again without loosing balance. This shows, that the simple cart-table model is a sufficient approximation for large quadrupeds performing locomotion tasks.

4) Avoiding kinematic limits: Attempting to cross the gap with a statically stable gait tends to overextend the legs, since large body motions are required to move the robot into statically stable positions. Dynamic motions allow us to keep the CoG closer to the center of all four feet, since stability can be achieved by appropriate accelerations, avoiding kinematic limits.

5) Stability despite irregular swing-leg sequences: Walking over the stepping stones demonstrates the ability of the controller to execute irregular swing-leg sequences in a dynamically stable manner (Fig. 8). Starting from a lateral sequence gait (LH-LF-RH-RF) the foothold sequence changes to traverse these irregularly placed stepping stones. Despite this, the produced CoG trajectory (colored solid line) is dynamically stable, since the ZMP (asterisk) is always kept inside the current support triangle. When comparing the actual (top) and desired (bottom) CoG trajectories, a tracking error is evident. By keeping the ZMP e.g. d = 6 cm away from the stability borders, we are robust even against these inaccuracies.

The whole body motion generator inserts four-leg-support phases (red section) only whenever it detects disjoint support triangles in the swing-leg sequence. While executing steps 1 and 2 (Fig. 8) no four-leg-support phase is necessary, because the triangles are not disjoint. Only after returning to swing the right-front leg, the robot requires a four-leg support phase for the ZMP to transition from the green (LH) to the yellow (RF) support triangle at (x, y) = (1.1, 0).

VII. CONCLUSION

We presented a dynamic, whole-body locomotion framework that executes footholds planned on-board. We showed, how a change in the environment causes the foothold generator to re-plan footholds on-line. We presented a whole body motion planner, which is able to generate a ZMP-stable body trajectory despite irregular swing-leg sequences to execute footholds dynamically. We showed how a combination of virtual model and floating-base inverse dynamics control can compliantly, yet accurately, track the desired wholebody motions. Real world experimental trials on challenging terrain demonstrate the capability of our framework.

We are currently working on bringing the kinematic planning and dynamic execution closer together. The idea is to produce desired state trajectories and required torques through one trajectory optimization problem, taking into account torque/joint limits, the dynamic model of the robot, foothold positions, friction coefficients and other constraints. With this approach we aim to produce even more dynamic motions such as jumping and rearing, during which fewer or no legs are in contact.

Fig. 1 .

 1 Fig.1. The hydraulically actuated and fully torque controlled quadruped robot HyQ[START_REF] Semini | Design of HyQ -a hydraulically and electrically actuated quadruped robot[END_REF]. The inset plot shows the on-line built environment perception alongside with the planned footholds and the current on-board robot state estimate.

Fig. 2 .

 2 Fig.2. An overview of the perception and planning system that generates footstep sequences according to the terrain information.

Fig. 3 .

 3 Fig.3. A sketch of the body action graph. The objective is to find a sequence of actions a from the current body state s = (x, y, θ) to the goal state g, that minimizes the accumulated action costs c(s, a). For simplicity only three possible actions are shown, namely move left (a l ), right (a r ) and forward (a f ). The optimal action sequence {a l , a r , . . . , r } found through ARA* is shown in red.

Fig. 4 .

 4 Fig. 4. (Re-)planning and perception. The left image shows how a map of the environment is built (cost values in grayscale) along with the body path (green line) and the footstep sequence plan (colored spheres). Once the environments changes the map is updated and the footsteps re-planned.

Fig. 5 .

 5 Fig. 5. Left: Cart-Tablemodelfor representing a quadruped robot: The total mass of the robot is concentrated in the cart that moves on the table. The base of the table represents the current area of support, determined by the current footholds. The ZMP must lie inside this area for dynamic stability. Right: Disjoint support triangles due to the added stability margin d. When switching between swinging the left-front (LF) to right-hind (RH) the ZMP must move from the brown to the blue support triangle. Since all four feet are in stance during this phase, the ZMP is only restricted by the red support polygon.

Fig. 6 .

 6 Fig.6. Pipeline that uses planned footholds to generate dynamic wholebody motions and compliantly executes them using a combination of feedforward and feedback terms.

Fig. 7 .

 7 Fig. 7. Snapshots of experimental trials used to evaluate the performance of our framework. From top to bottom: crossing a 15 cm pallet; climbing a stair-like structure consisting of two stacked pallets; traversing a 35 cm gap and crossing a sparse set of stepping stones.

Fig. 8 .

 8 Fig.8.Top: The body motion when walking over the stepping stones is show in black. The planned footholds are shown and the irregular step sequence LF(1) → LH(2) → RF(3) is highlighted (red). Bottom: The 3 shrunk support triangles corresponding to the highlighted step sequence brown → green → yellow are shown. Additionally the planned CoG (solid line) and ZMP trajectory for the duration of these 3 steps is illustrated (asterix). While the CoG (solid line) does not reach the support triangles, the ZMP does, causing dynamic stability. When switching between disjoint support triangles (green → yellow) four-leg support phases are inserted (red) to allow a smooth transition.

TABLE I FORWARD

 I SPEED AND SUCCESS RATE OF EXPERIMENTS AVERAGEDOVER 10 TRIALS AND COMPARED TO PREVIOUS RESULTS FROM[START_REF] Winkler | Path planning with force-based foothold adaptation and virtual model control for torque controlled quadruped robots[END_REF].

		Speed [cm/s]		Success Rate [%]
	Terrain	Curr. Prev. Ratio	Curr. Prev. Ratio
	Step. Stones	7.3	1.7	4.2	60	70	0.9
	Pallet	9.5	2.1	4.5	100	90	1.1
	Two Pallets	10.2	1.8	5.8	90	80	1.1
	Gap	12.7	-	-	90	0	-
	graph together with ARA* allows us to replan footholds at a
	frequency of approximately 0.5 Hz for goal states up to 5 m.
	2) Speed while dynamically stable: The pallet climbing
	and gap experiments show the speed (Table

http://youtu.be/MF-qxA_syZg

ACKNOWLEDGEMENTS

This research is funded by the Fondazione Istituto Italiano di Tecnologia. Alexander W. Winkler was partially supported by the Swiss National Science Foundation (SNF) through a Professorship Award to Jonas Buchli and the NCCR Robotics. The authors would like to thank the colleagues that collaborated for the success of this project: Roy Featherstone, Marco Frigerio, Marco Camurri, Bilal Rehman, Hamza Khan, Jake Goldsmith, Victor Barasuol, Jesus Ortiz, Stephane Bazeille and our team of technicians.


John Pachón-Morales

Julien Colin

Floran Pierre

Thibaut Champavert

François Puel

Patrick Perré

Flowability of lignocellulosic biomass powders: influence of torrefaction intensity

The poor flowability of powders produced from raw lignocellulosic biomass may be an economically issue for the production of second-generation biofuels. Torrefaction is a pre-treatment step of the gasification process that improves the physical characteristics of biomass by making it more coal-like. Particularly, the loss of resilience allows a reduction of the grinding energy consumption and is likely to improve the flow behaviour of woody powders. In this study, we investigated the effect of particle size and shape distribution on flow properties (unconfined yield stress and flowability factor) of powder from raw and torrefied biomass (Picea abies). Several intensities of torrefaction were tested, and its extent was quantified by the global mass loss, chosen as synthetic indicator of torrefaction intensity (its accounts for both the temperature level and the residence time). The intensity of torrefaction shifts the particle size distribution towards smaller sizes. An effect on the circularity and aspect ratio was also observed. A strong, positive correlation was obtained between the measured flowability of biomass powders at different consolidation stresses and the intensity of heat treatment. These results confirm the interest of torrefaction as a pre-treatment step and aim to provide new knowledge on rheological properties of biomass powders.

Introduction

In spite of the recent drop of oil price, the development of fuel from renewable resources is absolutely required to face global warming. After the first generation of biofuels that consumes edible products, a second generation using the lignocellulosic part of plantssuch as wood -, is yet to be developed at the industrial scale. One possible BtL (biomass to liquid) route includes a gasification process followed by a Fischer-Tropsch synthesis. Most of gasification processes require a supply of ground biomass in the form of powder for three main reasons: to ease the conveying of the raw material along the processing chain, to optimize its injection into the gasifier and to increase the chemical reactivity of the biomass. Currently, biomass transport is a key obstacle for the cost-competitive production of second-generation biofuels. Indeed, due to its fibrous structure, the grinding of native biomass is an energy intensive process and the powder obtained does not flow well [1]. Thus, a pre-treatment step, such as torrefaction, is needed to improve the process.

Torrefaction is a thermal treatment operating at moderate temperature (200 to 300°C) under atmospheric pressure in absence of oxygen. The process is traditionally characterised by low particle heating rates (<50 °C.min -1 ) and by a relatively long reactor residence time (typically higher than 30 minutes) [START_REF] Bergman | Petten[END_REF]. The resulting product is intermediate between wood and charcoal and exhibits several advantages when compared to the original material such as a decreased hygroscopicity, a dimensional stabilization, an increased resistance to biological decay, a loss of mechanical resilience, etc. [START_REF] Almeida | [END_REF] Several researchers have studied the question of estimating biomass flowability. Nevertheless, despite the frequent occurrence of handling problems, there is still a relatively poor knowledge of the flow properties of biomass materials [1]. More research is therefore needed in this field to ensure the sustainability of the supply chain [4]. This includes a better understanding of the flowability of lignocellulosic biomass powders, namely the combined effect of the torrefaction and grinding parameters. This is needed for engineers to develop cost-effective storing, conveying and processing strategies.

In this study, we investigate the effect of different torrefaction intensities on the flow properties of ground biomass, and assess the relations between size and shape of the particles and the flow behaviour of these powders.

Materials and methods

Material preparation

Wood Chip preparation

A spruce (Picea abies) tree from a plantation in the forest of Le Châtaignier in Riotord (France) was used in this work.

All samples came from the same board. The size of the samples was chosen to fulfil several constraints such as the amount required for flowability tests (width), the grinding capacity of the grinder (maximum length) and the homogeneity of the heat treatment (small thickness). The size of each individual wood chip was thus fixed to 20x20x3 mm 3 . Each sample number of Table 1 is composed of 32 chips of this size.

As it is well-known that the wood properties may vary in the tangential and radial directions, the 32 wood chips of each sample are selected following a recurrent stepped scheme along the board in order to average this variability and to obtain twin samples (Fig. 1).

Fig. 1.

Stepped cutting plan of the spruce board (Each colour corresponds to one individual sample of 32 chips). L: longitudinal direction, T: tangential direction.

Torrefaction procedure

A batch torrefaction device specifically developed in our laboratory was used to treat our sample [5]. The atmosphere of a hermetically closed Memmert UFP400 chamber was controlled by a nitrogen flow rate (5 L.min - 1 ) to reduce the oxygen content, thereby avoiding oxidation and ignition. In all cases, the oxygen level remained below 1.5%. A powerful fan inside the chamber insures an efficient heat transfer either to heat-up the sample or to limit thermal overshot due to exothermic reactions. The gas temperature was measured nearby the samples using a calibrated thermocouple.

A total of 6 treatments were conducted at 220, 250 and 280 °C for 1 or 5 hours according to the following protocol: a. heating from room temperature (about 20 °C) at a rate of 5 °C min -1 up to 100 °C. b. plateau at 100°C for 30 minutes in order to evaporate the residual bound water before the heat treatment step. c. heating at a rate of 5°C min -1 until the treatment temperature. d. plateau at the treatment temperature for the desired duration and e. cooling down thanks to thermal losses and the nitrogen flow entering the reactor.

After treatment, samples were weighed to determine the mass loss. Mass loss (ML) due to heat treatment was calculated as follows:

100 (%) 0 0    m m m ML t (1)
Where, m0 and mt are the oven-dry mass before and after torrefaction, respectively. Mass loss is known to be a good indicator of the treatment intensity [START_REF] Bergman | Petten[END_REF]. Non-treated 0

Grinding procedure

Heat-treated samples were ground in a M20 IKA semiautomatic laboratory 550W knife miller. The rotational speed of the knives was 20 000 rpm. The same grinding protocol was applied to all of the samples: a 95 s grinding made in several cycles [START_REF] Crawford | [END_REF] with intermediate waiting times of 15 s in order to avoid any additive heat treatment due to thermal dissipation. A circulation of cold water was applied to reduce heating during grinding.

Morphological granulometry

A Sympatec-QICPIC morphological granulometer with a M7 lens (measuring range from 10 to 10240 µm) with a GRADIS dispersion and VIBRI feeding units was used in this work [START_REF]Laser Diffraction Sensors[END_REF]. To conduct the test, a sample of about 5 g of wood powder was loaded onto the hopper of the instrument. The size and shape of powders were recorded and analysed by the manufacturer software. The Feret mean diameter, the mean value of the distance of two tangents to the particle contour over several orientations was measured: The mean diameter of the sample is then calculated as:

                 n i i n i i i p n d n d 1 1 / ) 0 , 1 ( (2)
Where ni is the number of particles with a Feret mean diameter di. Shape of the particles is characterised through the circularity and the aspect ratio. For the current study, circularity, ψS (from 0 to 1), is defined as the ratio of the perimeter of a circle having the same area as the projected area of a particle to its actual perimeter. A perfect sphere would give ψS =1 and the smaller the value of ψS, the more irregular is the shape of the particle. The aspect ratio, the ratio of the minimal to the maximal Feret diameter, gives an indication of the elongation of the particle.

Flow properties testing

A RST-XS Schulze ring shear tester was used to assess the flow properties of biomass powders. The ring shear tester is a widely used device to measure flow properties of powders, including angle of internal friction, wall friction and bulk density. The standard procedure leads to results with low variability [6].

A set of three pre-consolidation stresses (σpre), considered as representative of the stress range for industrial applications, were tested: 2, 5 and 10 kPa [START_REF] Schulze | Powders and bulk solids[END_REF]. Three shear points at 25, 50 and 75 % of σpre for each consolidation stress were used to determine the yield locus of each sample. Using Mohr circle analysis, the yield loci obtained allowed to evaluate the consolidation stress σ1 and the unconfined yield strength σc for the different normal stresses at preshear applied. The ratio FFC (flow function coefficient) of consolidation stress to unconfined stress (σ1/σc) is used to characterize flowability numerically. The larger FFC is, i.e., the smaller the ratio of the unconfined stress to the consolidation stress, the better a bulk solid flows.

Results and discussion

Particle Size distribution

In Fig. 2 the cumulative size distributions of the powder samples are represented. The intensity of heat treatment, represented by the mass loss of the samples, has an obvious impact on the size distribution of the biomass powders: the loss of resilience allows to obtain smaller particles by grinding.

Fig. 2. Cumulative passing at several torrefaction intensities.

The reduction of the resilience at increasing intensity of heat treatment makes the samples brittle, thus the particle size distribution shifts towards smaller particle sizes and its dispersity is reduced. Similar conclusions have been also observed from other studies, for different kinds of biomass powders and torrefaction conditions [10], [11], [START_REF] Arias | [END_REF].

Fig. 3 shows the value of the mean diameter and the cumulative percentage of particles with a Feret mean diameter below 101.2 µm at various torrefaction intensities. The cumulative number represents the fraction of particles having a diameter of 101.2 µm or less.

The mean diameter decreases dramatically whereas d101.2 increases with the increase of torrefaction intensity, the grinding conditions remaining the same for all the samples. This shows that there is an improvement in the grindability characteristics of treated biomass. This result is in good agreement with those in previous works [START_REF] Arias | [END_REF], [9]. Fig. 4 represents the morphological parameters of the biomass powders for two torrefaction intensities. On the whole, circularity of spruce powders was improved as the intensity of treatment increased from ML=11.4 % to ML=46.5 % for particles size above 50 µm. More spherical particles are thus globally obtained when a severe treatment is applied (280°C, 5h). This denotes a loss of the fibrous behaviour. However, no significant difference in circularity of ground biomass was observed for torrefaction intensities from 11.4 % to 30.9 %. It seems that the finer particles tend to be more spherical. Nevertheless, a more sensitive morphological analysis is necessary in order to better characterize these small particles (dp<100 µm). Similar conclusions can be drawn regarding the shape factor: more uniform particles are obtained when intensity of torrefaction is increased. dependence between the extent of the heat treatment and the flowability of spruce powders is depicted. The native ground biomass has a cohesive behaviour that decreases with the heat treatment until a nearly free flowing feature for the most intensively treated sample. In almost all cases, the unconfined yield stress increases with the consolidation stress to the same extent, so the regime of flowability for the powders is independent of the consolidation stress applied (except for the most intensely treated sample for which a free-flowing behaviour is reached at elevated consolidation stresses). For all samples, better flowability is obtained at a greater consolidation stress but there is no influence of the stress level on the qualitative ranking of the set of powders. This last feature is in good agreement with the general behaviour for most of the powders reported in literature [START_REF] Schulze | Powders and bulk solids[END_REF]. A direct relationship between the intensity of torrefaction and the flowability factor was found (Fig. 6). For sake of clarity, only the values at σpre= 10 kPa are presented, but similar tendencies were verified for the others values. The improvement of flow behaviour with the extent of heat treatment could be seen as a combined effect of several factors: the modifications on particle size and width of the distribution, the improvement of shape characteristics and the modification of the surface interactions. Further stage of our study will consist of isolating each factor in order to discriminate their effects on flowability. The lower values of circularity and larger aspect ratios lead to more frequent interparticle entanglements, which degrades flowability. These findings are in agreement with literature reports [START_REF] Crawford | [END_REF], [START_REF] Almendros | 19th European Biomass Conference[END_REF].

Flowability

Conclusions

This work evaluated the impact of six operating conditions of torrefaction (from 11.4 % to 46.5 % of mass loss) on the flowability of a biomass powder measured in a ring shear tester, as well as on size distribution and morphology. Several tendencies are clearly drawn from this work:

-the particle size distribution shifts towards smaller particle sizes and its dispersity is reduced with increase in torrefaction intensity, -due to the loss of the fibrous behaviour, the circularity of spruce powders is improved at high-treatment intensities, -an improvement of the flowability of spruce powders by heat treatment is clearly exhibited. This behaviour gradually shifts from a cohesive behaviour for the native biomass to a nearly free flowing behaviour for a heat treatment with 46.5 % of mass loss.

Fig. 3 .

 3 Fig. 3. Effect of intensity of heat treatment on the mean diameter and cumulative percentage of particles with a size below 101.2µm.

Fig. 4 .

 4 Fig. 4. Shape of the treated spruce powders at two different torrefaction intensities.

Fig. 5

 5 Fig. 5 presents the flow functions for different torrefaction intensities along with the flowability regimes from not flowing to free-flowing behaviour. A clear

Fig. 5 .

 5 Fig. 5. Instantaneous flow functions for several intensities of torrefaction and lines of constant flowability FFC (solid lines).

Fig. 6 .

 6 Fig. 6. Effect of torrefaction intensity on the flowability factor of the biomass powders. σpre = 10 kPa.

Table 1 .

 1 Heat treatment intensities.

	Sample number	Treatment temperature (°C)	Duration of treatment (h)	Mass loss (%)
	1	220	1	11,4
	2	220	5	14,7
	3	250	1	17,4
	4	250	5	25,4
	5	280	1	30,9
	6	280	5	46,5
	7