R Chauvin
email: remi.chauvin@cea.fr

L Bennani
email: lokman.bennani@onera.fr

P Trontin
email: pierre.trontin@onera.fr

P Villedieu
email: philippe.villedieu@onera.fr

An Implicit Time Marching Galerkin Method for the Simulation of Icing Phenomena with a Triple Layer Model

Keywords: icing, triple layer, finite element method, runback, modelling, numerical simulation

In the context of more electrical aircraft and reduction of fuel consumption, aircraft manufacturers are moving towards more complex and transient ice protection systems. The operating of these systems involves several unsteady heat and mass transfer phenomena. Modelling and numerical simulation play an important role in the investigation of these unsteady phenomena. In this paper, a model for unsteady ice build-up and melting is presented. The model is based on a triple layer assumption. In addition, a tailored numerical methodology for solving the governing partial differential equations is also described. It is based on a Galerkin finite element method and a Gauss-Seidel like implicit time marching scheme. The global method is validated and its capabilities are demonstrated on several cases.

Introduction

Ever since the pioneering days of flight, the icing phenomenon, or the "ice problem" as it was called in the early days, has been reckognized as a serious threat [30,31]. Typical icing conditions are due to the presence in clouds of supercooled water droplets. Upon impacting an aircrafts surface, the metastable supercooled state of the droplets is broken. At that point, the droplets undergo a liquid-solid phase change, leading to ice bluid-up on the impinged surface. Increased mass, degradation of aerodynamic performances and handling qualities or blocked air intakes, those are some of the many undesirable consequences one has to face when dealing with icing. Indeed, in the aeronautical world, icing is one of the most serious hazards that can be encountered.

Aircraft manufacturers must therefore comply with certifications and regulations regarding flight safety in icing conditions. In order to achieve that goal, several ice protection technologies may be used. A commonly employed system for large airliners is the bleed-air anti-icing system. However, in the context of more electrical aircraft and reduction of fuel consumption, alternative systems are currently under consideration for large airliners. The main idea is to use a system which would be able to function in de-icing mode. That is to say that a reasonable amount of ice is allowed to build up. Periodically, the system is activated and the ice is shed from the surface or melted. Such a system would have an unsteady cyclic operating mode, enabling a more optimized use of energy. These systems operate in complex environments and involves many physical processes. In order to study such systems and assist their design phase, models and numerical simulation methodologies are an important asset.

Historically, one of the first models devised to study ice build up on a surface was Messinger's model [17]. It has served as the basis for many icing simulation tools [12,[START_REF] Gent | TRAJICE2-A combined water droplet trajectory and ice accretion prediction program for aerofoils[END_REF][START_REF] Dillingh | Simulation of ice accretion on airfoils during flight[END_REF][START_REF] Hedde | ONERA three-dimensional icing model[END_REF][START_REF] Paraschivoiu | Ice accretion simulation code canice[END_REF]32] and is still in wide use to this day.

However, it is a steady state model and is not well suited if one wishes to take into account the interaction with an unsteady system. Moreover, an important underlying hypothesis to the model is that the temperature does not vary in the ice block in the direction normal to the surface. As a consequence the heat fluxes are not taken into account and the model is not well suited to applications involving ice protection systems. For example, Messinger's model is not able to capture the melting dynamics of an ice block when an ice protection system is activated. In addition, the operating of ice protection systems can involve freezing runback films which are not well accounted for by Messinger's model.

In order to model and simulate unsteady icing phenomena, several issues must be addressed. First, one must account for the dynamics of the liquid film. This issue has been studied by many researchers. For example, Al-Khalil et al. have proposed a macroscopic runback model applied to the study anti-icing systems [START_REF] Km Al-Khalil | Further development of an anti-icing runback model[END_REF]. In addition, Tsao and Rothmayer have performed theoretical work using matched asymptotic expansions [28,27]. In particular, their analysis gives a possible explanation to the formation of roughness at the surface of glaze ice. Furthermore, in order to perform unsteady simulations, Bourgault et al. proposed an unsteady three dimensional ice accretion model based on a lubrication theory assumption [START_REF] Bourgault | Development of a shallow-water icing model in fensap-ice[END_REF]. Although more sophisticated approaches have been developed and researched [21, 10], lubrication theory is still widely used and sufficient for most applications [29,[START_REF] Lavoie | Numerical algorithms for infinite swept wing ice accretion[END_REF].

Secondly, one must take into account the unsteady heat and mass transfer phenomena. To do so, Myers has proposed to extend Messinger's model by using a multi-layered approach (one layer of ice and two layers of liquid water) [18,20]. In each layer, a linear temperature profile is assumed (in the direction normal to the surface) and the evolution of the temperature field is treated in a quasi-steady framework. These ideas were used in more recent work in a quasi-steady framework [29,33,11]. Later, Myers et al. extended this idea to unsteady one dimensional cases by assuming a cubic temperature profile in each layer [19]. The cases addressed in [19] were of fixed total thickness. Phenomena such as evaporation or ice build up due to impacting droplets were not treated. Moreover, the problem of creation/vanishing of a layer was not addressed.

The goal of the present work is to extend the previously mentioned research to develop a unified model and numerical method, enabling to take into account all unsteady phenomena encountered during ice buildup on a surface (eventually heated by a system). The proposed approach allows to take into account the necessary unsteady physics (water runback, dynamic phase change interfaces, unsteady temperature fields, creation/vanishing of a layer, etc) while remaining computationally tractable. Moreover, the method allows numerical coupling with simulation tools for thermal ice protection systems.

The present paper is articulated as follows. In the first section, the general context of the approach is presented. After that, the main physical models, which constitute the building blocks of the three-layer approach, are introduced. Then, the spatial finite element discretization and implicit time marching scheme are described. Finally, some validations and applications of the method are shown.

General Context

The simulation of icing phenomena in aeronautics involves several coupled physical phenomena. State of the art icing codes are articulated around specialized modules which are coupled during the computation.

For example, a widely used approach to compute an ice shape is illustrated in Figure 1. An inviscid flow solver is used to obtain the external flow field. This field is used as an input for a droplet trajectory module and a boundary layer module. These two modules in turn provide the inputs for the ice accretion module.

Inviscid flow solver

Boundary layer solver

Droplet trajectory solver

Ice accretion solver

MiLeS

Multi-Layer Ice Solver -Running film -Ice shape -Melted areas As stated earlier, classic ice accretion solvers are based on Messinger's model. The idea is to partition the airfoil into cells and to consider each cell as an open thermodynamic system. The equations of conservation of mass and energy in steady state are solved, but by considering a uniform state on the cell. As shall be described in the rest of the paper, the model introduced here is unsteady and considers three layer so as to account for the running film dynamics and the temperature gradients within the ice and the melted water regions.

The work presented here is focused on the 'ice accretion solver' block of Figure 1. The unsteady ice accretion solver is called MiLeS and its numerical implementation is part of a larger simulation framework, developed at ONERA, called IGLOO2D and IGLOO3D [22,24]. In what follows, the model and numerical method will be presented in a two dimensional framework.

As shown in Figure 1, the inputs of the ice accretion solver are the (local) droplet impingement rate noted ṁimp , their total temperature T tot imp , the heat transfer coefficient h tc , the recovery temperature T rec , the shear stress τ and the inviscid flow quantities. The ice accretion solver outputs the local running film height, height of ice, height of melted water and the temperature fields in each layer.

Physical Model of MiLeS2D

In order to construct the three-layer model, let us consider the illustrative icing situation shown in Figure 2. In this situation, the physical process may be described by considering six distinct modes (labeled 1 -6 in Figure 2):

1. Full evaporative: the whole mass of impacting droplets is evaporated (for example, due to heat provided by an ice protection system). 2. Running wet: only a liquid water film is present. Under the action of the aerodynamic forces, the liquid film runs back along the surface. 3. Rime accretion: the droplets freeze almost instantaneously leading to ice build up with no liquid water. 4. Glaze accretion: the droplets freeze, but at a slower rate than in the rime case. Therefore, a running liquid water film is present on top of the ice layer. 5. Rime accretion with melting at the surface (due to heat provided by an ice protection system for example). 6. Glaze accretion with melting at the surface (due to heat provided by an ice protection system for example).

Ice

Liquid water 1 In each mode, one can distinguish the presence (or abscence) of one or more layers. The term three-layer model comes from the fact that in the most general case one has to deal with three layers:

1. Running film: this is the water layer directly in contact with the airflow. It is sheared by the aerodynamic flow and runs back downstream. 2. Ice: this layer represents the ice that may build up in icing conditions. 3. Static film: when a system is activated, a part of the ice layer may melt. This has the effect of forming a static film layer (it is not in contact with the external aerodynamic flow).

The three layers, as well as some geometric quantities which characterize them, are represented in Figure 3. Subscript f is used to denote the running film layer, i for the ice layer and s for the static film layer. For each layer k (k = i, s), a k is the position of the lower boundary and b k is the position of the upper boundary (in the z direction, normal to the surface). The whole model may be fully expressed in terms of PDEs [7]. However, this requires to consider each mode shown in Figure 2 because specific boundary conditions and source terms correspond to each one. In order to present the approach in a more concise way, while not losing much in generality, it shall be presented as follows. First, the main equations used to model heat and mass transfer will be presented. Secondly, the ways of combining them to describe each of the six modes are given with the corresponding set of boundary conditions and source terms. Then, in the following section, the spatial discretization will be presented. Finally, the implicit time integration scheme is described.

Static Film

Equations for the running film layer

In order to obtain an unsteady formulation, the dynamics of the liquid film need to be modelled. The unknowns are the film thickness h f and the mean temperature Tf . To do so, as the film is very thin, an integral approach is used combined with a closure relation from lubrication theory [7].

In the framework of lubrication theory, the mean velocity of the running liquid (v x ) only depends on its thickness h f and the external flow conditions:

vx (h f ) = τ 2µ w h f + 1 3µ w - ∂p ∂x + ρ w g x h 2 f ( 1
)
where τ is the air shear stress, ∂p ∂x is the tangential pressure gradient and g x is the tangential gravity. In this framework, surface tension and wettability effects are neglected.

One obtains the following system of PDEs:

• Mass conservation

∂ρ w h f ∂t + ∂ρ w h f vx ∂x = Γ f (2)
where Γ f is the mass transfer term (which will be detailed later and is specific to each mode).

• Energy conservation: as the liquid film is very thin, the temperature profile is assumed uniform in the normal direction. The energy equation is given by:

∂ρ w c w h f Tf ∂t + ∂ρ w c w h f vx Tf ∂x = Φf ( 3
)
where Φf is the energy transfer term (which will be detailed later and is specific to each mode).

In the case where the liquid film is (locally) running over ice, equation ( 3) degenerates to Tf = T m . Indeed, at an interface between running film and ice, phase change is occurring. Hence the temperature at the interface is the melting temperature. As the liquid film is assumed to be very thin and of uniform temperature on the normal direction, the whole film is at the melting temperature.

Heat conduction in the ice layer and melted film layer

Both the ice and melted film layers are assumed motionless. Here, the subscript k will be used to generically denote these two layers (k = i or k = s). Furthermore, heat transfer in the tangential direction is neglected (see Appendix A for an assessment of this approximation). Heat transfer is considered in the z (normal) direction only.

• Mass conservation (at each curvilinear abscissa x):

∂ρ k h k ∂t = Γ k (4)
where Γ k is the mass transfer term (which will be detailed later).

• Energy conservation (at each curvilinear abscissa x)

∂ρ k c k T k ∂t = λ k ∂ 2 T k ∂z 2 For z in ]a k (t), b k (t)[ f a T k , ∂T k ∂z = 0 For z = a k (t) (5)
f b T k , ∂T k ∂z = 0 For z = b k (t)
where for each layer k = i or s, ρ k is the density, c k the specific heat and T k the temperature. f a and f b denote in a generic way the boundary conditions at the moving boundaries of each layer, as explained in the next paragraph.

The six modes: boundary conditions, source terms and coupling between layers

In the previous subsection, the main governing equations for each layer (running film, ice or static film) were described. In order to obtain a complete model, these equations are combined so as to describe the six possible modes shown in Figure 2.

Before describing each mode, let us first discuss some general aspects. The convective heat transfer with the outer boundary layer and internal surface heating are independent of any given mode. The convective heat transfer flux is given by the relation:

Φconv (T ) = h tc [T rec -T ] ( 6
)
where h tc is the boundary layer heat transfer coefficient and T rec is the recovery temperature. They are inputs of the model.

The heating flux at the lower boundary is given by

Φwall (T ) = h wall [T int -T ] + Φ0 (7)
where h wall is the wall heat transfer coefficient, T int is the temperature on the inside of the surface and Φ0 a given flux (in the case of a Neumann condition for example). These are again inputs of the model. In what follows, ṁev denotes the evaporation mass rate and ṁs the sublimation mass rate. They are both computed using the Chilton-Colburn analogy [START_REF] Chilton | Mass transfer (absorption) coefficients prediction from data on heat transfer and fluid friction[END_REF].

At interfaces where phase change occurs the rate of phase change is governed by a Stefan condition. For example, at the interface between the ice and the static water layers, ṁice (rate of phase change of the ice layer) is given by:

ṁice L m = λ w ∂T s ∂z interf ace -λ i ∂T i ∂z interf ace (8)
In addition, recall that the temperature profile in the running film is assumed to be uniform. In the modes were a liquid running film layer is present, the computation of the temperature gradient in the running liquid film layer is not possible. However, this term is required in order to compute the rate of melting/freezing (governed by the Stefan condition as just stated). So as to overcome this limitation, the temperature gradient is also assumed constant (and different from zero) as if the temperature profile was in fact linear.

Note that the uniform temperature and constant gradient assumptions are not incompatible in the present case. The uniform temperature assumption is a closure hypothesis for the energy equation for the running film [7]. It is a reasonable assumption given the small characteristic thickness of the running film in the targeted applications of this work. The characteristic thickness of the running film is of a few microns. A temperature difference of 0.1K across the film thickness, which is coherent with a uniform temperature assumption, would yield a temperature gradient of 10 5 K/m (assuming a linear profile). Therefore, the uniform temperature and constant gradient assumptions are coherent (in our context).

When the heat flux at the running film/ice interface is required, it is computed by assuming that it is equal to the heat flux at the upper boundary of the running film. This is used for example in equation ( 17) for the glaze mode.

We now proceed to the description of each mode and their specific boundary conditions, source terms and coupling terms.

-Full evaporative

In this mode, the whole quantity of water in the cell is evaporated. No layer is present in this mode. The unknown is the temperature at the wall T f e . The full evaporative case is governed by equations 3) are given by:

Γ f = ṁimp -ṁev ( Tf ) (11) Φf = ṁimp c w T tot imp + ṁev ( Tf ) c w Tf -L v ( Tf ) + Φconv ( Tf ) + Φwall ( Tf ) (12)
3 -Rime accretion This mode corresponds to the case where all the water arriving in the cell freezes instantaneously. The unknowns are the ice layer height h i and the temperature field in the ice T i . The equations describing this mode are ( 4) and ( 5). The source terms and boundary conditions are:

Γ i = ṁimp -ṁs (T i (h i )) - ∂ρ w h f vx ∂x ( 13
)
-λ i ∂T i ∂z z=0 = h wall [T wall -T i (0)] + Φ0 (14)
λ i ∂T i ∂z z=hi =h tc [T rec -T i (h i )] -ṁs (T i )L s (T i ) + ṁimp c w T tot d -T m -c i [T i (h i ) -T m ] + L f (15) - ∂ρ w c w h f vx ( Tf -T m ) ∂x - ∂ρ w h f vx ∂x [c i [T i (h i ) -T m ] + L f ]
where the two last terms of equation (15) represent respectively the energy flux brought by the running film and the energy accounting for the freezing of this mass of water.

-Glaze accretion

This mode describes the case where an ice layer is present with a liquid running film on top. The unknowns are h f , h i , and T i . As in this mode, there is a liquid layer on top of the ice layer, equation (3) degenerates to Tf = T m . This mode is governed by equations (2), ( 4) and ( 5) with the following source terms and boundary conditions:

Γ f = ṁimp -ṁev (T m ) -ṁice,top (16)
where ṁice,top is the rate of phase change at the interface between the running liquid film and the ice.

ṁice,top

L m = λ i ∂T i ∂z z=hi -ṁimp c w T tot imp -T m + Φconv -ṁev (T m )L v (T m ) - ∂ρ w c w h f vx Tf ∂x (17)
where the last four terms are due to the uniform temperature gradient assumption in the running film flowing over a block of ice.

Γ i = ṁice,top (18)
T i (h i ) = T m (19) -λ i ∂T i ∂z z=0 = h wall [T wall -T i (0)] + Φ0 (20)

-Rime accretion with static film

In this mode, two static layers are present, the ice layer and a static melted film layer. The unknowns are h i , T i , h s and T s . This mode is governed by equations ( 4) and ( 5) with the following source terms and boundary conditions:

Γ i = ṁimp -ṁs (T i (h s + h i )) + ṁice,bottom - ∂ρ w h f vx ∂x (21)
where ṁice,bottom is the rate of phase change at the interface between the static layer and the ice.

ṁice,bottom

L m = λ w ∂T s ∂z z=hs -λ i ∂T i ∂z z=hs (22)
T i (h s ) = T m (23)
λ i ∂T i ∂z z=hs+hi =h tc [T rec -T i (h i )] -ṁs (T i )L s (T i ) + ṁimp c w T tot d -T m -c i [T i (h i ) -T m ] + L f (24) - ∂ρ w c w h f vx ( Tf -T m ) ∂x - ∂ρ w h f vx ∂x [c i [T i (h i ) -T m ] + L f ] Γ s = -ṁice,bottom (25)
-

λ s ∂T s ∂z z=0 = h wall [T wall -θ s,1 ] + Φ0 (26)
T s (h s ) = T m (27)
6 -Glaze accretion with static film This is the three-layer mode, a film is running on top of an ice layer with a static melted film layer underneath. The unknowns are h f , h i , T i , h s and T s . This mode is governed by equations (2), ( 4) and ( 5) with the following source terms and boundary conditions:

Γ f = ṁimp -ṁev ( Tf ) -ṁice,top (28)
Γ i = ṁice,top -ṁice,bottom (29)
T i (h s + h i ) = T m (30) T i (h s ) = T m (31) Γ s = -ṁice,bottom (32)
-

λ s ∂T s ∂z z=0 = h wall [T wall -T s (0)] + Φ0 (33)
T s (h s ) = T m (34)

Numerical Method

The next step is to discretize the problem at hand so as to numerically obtain a solution. The spatial discretization is performed as follows:

i The airfoil is discretized into N x control volumes. ii Using this x-wise mesh, the spatial derivatives in equations ( 2) and ( 3) are discretized with a finite volume Roe scheme [START_REF] Roe | Approximate riemann solvers, parameter vectors, and difference schemes[END_REF]. Moreover, in order to obtain a scheme enforcing the maximum principle, a method inspired by the one proposed by Larrouturou is used [15]. These schemes are well established and will not be described in further detail. iii The z-wise derivatives in the (eventual) ice layer or static film layer are discretized using a Galerkin method.

In the following subsections, the spatial discretization for the heat equation and the implicit time integration scheme are described.

Finite element spatial discretization in the z direction for the heat equation

If in a given mode, the cell contains an ice layer and/or a static film layer, the heat equation [START_REF] Bennani | Numerical simulation and modeling of ice shedding: Process initiation[END_REF] needs to be solved. The main difficulty here is that this equation is written in a moving interface domain delimited by a k (t) and b k (t). The change in boundary position is due to phase change phenomena, impingement and evaporation.

The first step is to write the variational formulation of (5). To do so, let v(z) be a test function. The variational formulation of (5) reads:

b k a k ρ k c k ∂T k ∂t v dz - b k a k λ k ∂ 2 T k ∂z 2 v dz = 0 (35)
The second step of this method is to introduce the change of variables:

z(t, z) = 2(z -a k (t)) b k (t) -a k (t) -1 (36)
This allows to define the problem on [-1, 1] instead of [a k , b k ]. Under this change of variables, the variational formulation takes the form:

1 -1 ρ k c k ∂T k ∂t + ∂ z ∂t ∂T k ∂ z v b k -a k 2 dz - 1 -1 λ k ∂ 2 T k ∂ z2 v 2 b k -a k dz = 0 (37)
The third step is the discretization of (35) in a subspace of finite dimension. Let E n = (Ψ j ) n j=1 be a set a linearly independent functions. The temperature field is projected onto this subspace:

∀t ≥ 0, ∀z ∈] -1, 1[, T k (t, z) = n j=1 θ k,j (t)Ψ j (z) (38)
Then one can project the equation ( 35) onto E n . This yields, ∀i = 1, . . . , n:

ρ k c k n j=1 1 -1 θk,j Ψ j (z) + ∂ z ∂t θ k,j Ψ j (z) Ψ i (z) b k -a k 2 dz -2λ k b k -a k n j=1 1 -1 ∂ 2 ∂ z2 [θ k,j Ψ j (z)] Ψ i (z)dz = 0 (39)
By developing this expression and using integration by parts for the second term, one gets:

ρ k c k n j=1 θk,j b k -a k 2 1 -1 Ψ j (z)Ψ i (z) dz + ρ k c k n j=1 θ k,j 1 -1 Ψ j (z) ∂ z ∂t Ψ i (z) b k -a k 2 dz -   2λ k b k -a k n j=1 θ k,j Ψ j (z)Ψ i (z)   1 z=-1 + 2λ k b k -a k n j=1 θ k,j 1 -1 Ψ j (z)Ψ i (z) dz = 0 (40)
By using this procedure one obtains the following matrix expression (details can be found in [7, p.78]):

ρ k c k b k -a k 2 A θk -1 2 ( ȧk + ḃk )Π + ( ḃk -ȧk )Γ θ k + 2λ k b k -a k Λθ k = f k (θ k ) (41)
where

                                             θ k = (θ k,1 . . . , θ k,n ) T A ij = 1 -1 Ψ j (z)Ψ i (z) dz Λ ij = 1 -1 Ψ j (z)Ψ i (z) dz Π ij = 1 -1 Ψ j (z)Ψ i (z) dz Γ ij = 1 -1 zΨ j (z)Ψ i (z) dz f k,i (θ k ) = 2λ k b k -a k n j=1 θ k,j Ψ j (1)Ψ i (1) - 2λ k b k -a k n j=1 θ k,j Ψ j (-1)Ψ i (-1) (42)
Note that the matrices A, Π, Γ and Λ only depend on the functions Ψ j . So as to reduce the computational cost they are evaluated once at the beginning of the computation. Moreover, f k,i (θ k ) is computed using the boundary conditions.

Finally, the functions Ψ j remain to be chosen. In this work polynomial functions are considered. They are chosen so that the first and second degree of freedom yield the value of the temperature at the lower and upper boundaries respectively. In addition they are constructed so as to be orthonormal. The constructions of those polynomials reads:

1. ∀x ∈ [-1, 1], Ψ 1 (x) = 1 -x 2 , Ψ 2 (x) = 1 + x 2 2. For i ≥ 3, Ψ i is the polynomial of degree i -1 such that (a) Ψ i (-1) = Ψ i (+1) = 0 (b) 1 -1 (Ψ i (z)) 2 dz = 1 (c) 1 -1 Ψ i (z)Ψ j (z)dz = 0, ∀j = 3, . . . , i -1
One may note that conditions 1. and 2.(a) allow to easily manage a Dirichlet boundary condition since

T k (a k ) = θ k,1 and T k (b k ) = θ k,2 .

Implicit time marching algorithm

To insure robustness and a lower computational cost, an implicit Euler time marching scheme is used. The time interval is discretized into N t discrete instants and superscript n is used to denote the n-th time step.

Let's assume that the computational domain has been split into N x cells (as illustrated in Figure 4). Then, for a given cell k in mode md, the general shape of equation to solve to go from time step n to time step n + 1 reads:

M md U n+1 k -U n k ∆t + G md (U n+1 L , U n+1 k , U n+1 R ) = S md (U n+1 k ) ( 43
)
where for all k ∈ 1, N x , the unknown vector is denoted U n+1 k . As an example, in the glaze mode, U n+1 k = (h f , h i , θ i,1 , ..., θ i,ng ) (where ng is the number of degrees of freedom of the Galerkin method). Subscripts L and R denote the left and the right neighbouring cells of cell k, respectively. G md is a nonlinear function defined by the spatial discretization of the running liquid film (Roe finite volume scheme). Finally, S md is a local function (it only depends on the variable inside the local cell) which models heat transfer in the normal direction including phase change as well as source terms (previously described Galerkin method). It is worth noting that the dimension of U n+1 k depends on k, md and n since layers may appear or disappear during a time step.

There are two main difficulties in the construction of the implicit method. First, as all the cells of the mesh are coupled with their left and right neighbouring cells through the running liquid film, it is complicated to devise a global implicitation of the method (i.e. a method which would simultaneously yield U n+1 k in every cell). To overcome this problem, a fixed point algorithm is used. The general idea is to perform a local implicitation and iterate over the cells until convergence so as to construct a sequence U (l) k k∈ 1,Nx l that converges to U n+1 k k∈ 1,Nx . More precisely, the first step is to chose a numbering for iterating over the cells. In order to enhance convergence, the idea is to run through the cells in the direction of propagation of information. It is given by the direction of shearing of the film, which determines the main direction of propagation of the film downstream. For the present targeted icing applications, the numbering starts at the stagnation point (in red in Figure 4). Then one firstly numbers the lower surface and then the upper one, as illustrated in Figure 4. The left and right neighbouring cells relative to a given cell k are defined with respect to this numbering direction.

M md U (l+1) k -U n k ∆t + G md U (l+1) L , U (l+1) k , U (l) R = S md U (l+1) k (46) If (U (l+1) k
is not a compatible mode (as defined by Figure 5)) then Change mode according to Figure 5 and solve (46) again end If until (Compatible mode is found)

end For until (∀k ∈ 1, N x , U (l+1) k -U (l) k U (l+1) k < ε)
Algorithm 1: Fixed point for the implicitation of time integration

Validation: The Stefan Problem

In this section, a Stefan problem [2] is considered so as to assess and validate the Galerkin method. Consider a block of ice, initially at 250K, as in Figure 6. The block is heated for 5 s from beneath with an imposed temperature of T L = 300 K. As the imposed temperature is greater than the melting temperature (273.15K) the block of ice will undergo a solid-liquid phase change and a melting front will propagate from beneath.

Liquid Solid

Figure 6: Illustration of the one dimensional Stefan problem: a block of ice is heated from beneath at a temperature T L . Both lateral sides are adiabatic (q = 0).

For a semi-infinite block, this problem has an analytical solution [2]. However, when performing the numerical simulation, the ice block has a finite size. Therefore, the length of the block is taken large enough so that the main heat transfer is confined to the bottom of the domain. In the present case the length of the block is 0.01m. The physical properties used for this case are given in table 1. The previously described method is used to solve the Stefan problem. Therefore, the initial mode is "rime". When the melting front starts to propagate, the mode switches to "rime with static film". The problem is solved using 2, 4, 6 and 8 degrees of freedom in the Galerkin expansion and a time step of 0.01s. Figure 7a shows the melting front position as a function of time. With 2 degrees of freedom, the numerical solution is overestimated. The solution is improved with increasing number of degrees of freedom. With 4 or more degrees of freedom the results are in good agreement with the analytical solution.

c s (J.kg -1 .K -1 ) c l (J.kg -1 .K -1 ) L f (J.kg -1 ) λ s (W.m -1 .K -1 ) λ l (W.m -1 .K -1 ) T m (K)
In addition, Figure 7b shows the sensitivity of the solution obtained with 8 degrees of freedom to the time step. The agreement with the analytical solution is very good. The solution obtained with ∆t = 0.1s is slightly less accurate. Solutions obtained with ∆t = 0.01s and ∆t = 0.001s are almost identical. So as to get a better idea of the validity of the method, one may also compare the numerical and analytical temperature fields. Here, the temperature fields are compared at the very first instants of the melting front propagation (t = 0.1s) and at a later time (t = 5s). As shown in Figure 8b, at t = 5s, except when using only 2 degrees of freedom, the agreement between analytical and numerical solutions is excellent. However, Figure 8a shows that at t = 0.1s, the numerical temperature fields tend to oscillate. This is due to the truncation operated by the Galerkin procedure. These oscillations tend to be lessened when a higher number of degrees of freedom is used. With 8 degrees of freedom, the results are in good agreement with the analytical solution. In the rest of this work, 8 degrees of freedom will be used for the Galerkin method.

Ice Accretion Cases

This section describes the results obtained for two steady ice accretion cases (rime and glaze). They are compared with ice shapes predicted with Messinger's method (implemented in the module MESSINGER2D of the ONERA icing tool IGLOO2D [22], the present method is implemented in MiLeS2D). In both cases, the airfoil is discretized using 128 points.

The icing conditions for both cases are given in Table 2. Rime icing conditions correspond to the case where the external temperature is well below 273.15K. In this case the supercooled droplets freeze almost instantly upon impacting the surface. Therefore in rime conditions, there is no water running back. In glaze icing conditions, the external temperature is below but closer to 273.15K. In this case, the droplets start to freeze upon impact but a certain amount of water runs back downstream.

Case

Airfoil Chord α ( Figure 9 shows that the ice shapes predicted with MiLeS2D and MESSINGER2D match very well. This section shows that in steady icing conditions, the present numerical methodology is able to predict the same results as Messinger's method. This section therefore provides a successful benchmark (for steady configurations) for MiLeS2D before performing a case for which the method is really designed (unsteady and heated leading edge).

• ) P ∞ (P a) T ∞ (K) M ∞ rH ∞ T W C (g.m -3 ) M V D (

Delayed Activation With Runback

We now proceed to show that the method is also able to treat unsteady situations of high interest for icing applications, which can't be treated with Messinger's approach.

Consider the glaze icing conditions of table 3 and assume that after 20s of accretion, an ice protection system is activated. The heat provided by the system is represented by an internal heat transfer coefficient of 1500W.K -1 .m -2 and an internal temperature of 323.15K. The system is active between s = -0.015m and s = 0.015m. This corresponds to an illustrative case of delayed activation. Some ice has already had the time to build up and an internal heating will be activated to eliminate it. Such a case can't be simulated using Messinger's approach. Indeed, it can't account for the dynamics of the melting process nor can it account for the dynamics of the running back water film created when all the ice is melted.

Airfoil

Chord During the first 20s of the case the internal heating is not activated. Unsteady phenomena occur during this step but only at the very beginning of the build up process.

α ( • ) P ∞ (P a) T ∞ (K) M ∞ rH ∞ T W C (g.m -3 ) M V D (
Figure 10 shows, at t = 0.25s, the liquid running film and ice thickness as well as the droplet impact rate as functions of the curvilinear abscissa. Different areas of accretion mode (rime or glaze) are visible. Between approximately s = -0.01m and s = 0.0m droplets are impacting, ice is building up and a liquid film is flowing on top (glaze mode). On the other hand, for s < -0.02m or in the neighborhood of s = 0.0m, ice is building-up in rime mode. : Droplet impact rate, running film and ice build up modes at t = 0.25s. The plotted thicknesses dimensionless. The maximum thicknesses are h rf,max = 1.9 µm for the running film and h i,max = 6.9 µm for the ice.

As shown in Figure 11 the film flows downstream under the effect of the aerodynamic forces. The film flows into regions where ice was growing in rime mode. Therefore the build up process changes downstream as the liquid film runs back and flows over the ice. Regions which initially started in rime mode switch to glaze mode. In addition, note that after t = 1s steady state is reached regarding the runback of the liquid film (whereas ice continues to grow). This unsteady process occurs only during the first instants of the build up process. As shown in Figure 12, this leads to build up of ice at the leading edge. After 20s the maximal ice thickness is of 0.84mm. From this point on a heating term is activated in the red dotted region. Note that for representation reasons the dotted region is not plotted exactly at the boundary. However, it is indeed at the upper boundary of the airfoil that the heating term in applied. That is to say, the heating term is directly applied as a (lower) boundary condition of the triple layer model. 7.2.

Step 2: activation of heating At t = 20s the internal heating is activated. The layer of ice starts to melt. Figure 13 shows the temperature field in the ice block at the stagnation point just after the internal heating is applied. The melting front propagates and a temperature gradient is clearly visible.

Airfoil Airfoil Airfoil

Air

Air Air Figure 14a shows the temperature profile and the position of the melting front along the normal direction to the stagnation point. The temperature gradient is clearly visible. The ice melting process is also visible in Figure 14b. The ice thickness decreases as it melts and is replaced by a layer of static film. Once the ice layer has completely melted the static layer comes into contact with the running film layer. At this point the mode switches to "running wet" and the whole thickness of the static layer is transferred to the dynamic layer. The dynamic layer is then very rapidly evacuated downstream and returns to its steady thickness. In addition, during the unsteady melting process, small region of melted water are released and flow downstream, as shown in Figure 15. The whole process occurs as follows:

• t = 20.03s: in a given cell at the end of the ice block, once the melting front has propagated through the whole thickness of ice, the local mode switches to 'running wet'.

• t = 20.08s: The liquid water comes into contact with the external flow field. The liquid film runs back downstream under the effect of the aerodynamic forces.

• t = 20.21s: however, as the liquid leaves the protected area, it is no longer heated enough. It starts to freeze under the effect of the cold external conditions.

• t = 20.3s: therefore, a small ice ridge starts build up outside of the heated area.

Static film grows as ice block continues to melt

A i r f o i l A i r f o i l A i r f o i l A i r f o i l

Step 3: formation of ice ridges

After 10s of internal heating, the protected area is free from ice. The internal heating manages to keep the leading edge free of ice. However it is not sufficient to be fully evaporative. Figure 16 shows the predicted ice shapes after 50s. The temperature field in the running liquid film is shown. The water is heated along the protected area. However, once the water exits the protected area, the cold external conditions rapidly lead to freezing. Therefore the water runs back and freezes further downstream, generating two runback ice ridges, one on the pressure side and the other on the suction side.

Conclusion and Perspectives

In this paper, a three-layer model for the simulation of unsteady icing was proposed. The methodology models the dynamics of the sheared liquid film with a lubrication assumption. The heat equation is used to model heat transfer in the normal direction in the static layers. Moreover the dynamics of the melting front is taken into account with a Stefan condition. The heat and mass transfer equations for the running liquid film are solved using a Roe scheme. In the static layers, a Galerkin finite element method is used to solve the heat equation and treat the Stefan condition. These spatial discretization schemes are embedded into an implicit time marching scheme.

The approach was shown to give good results when compared to analytical solutions for a one dimensional Stefan problem. It was shown that the melting front and temperature fields could be well predicted.

Secondly, two ice accretion cases where simulated. The results were compared to those obtained with Messinger's method and to experimental data. In both cases, rime and glaze icing conditions, the present method predicted the same results as Messinger's approach. However, results were also in poor agreement with the experimental data. In order to overcome this drawback, a possible line of improvement would be to update the aerodynamic field and droplet impingement data. As the present method is fully unsteady, one could update the aerodynamic field and droplet impingement data in a consistent way.

Finally, the capabilities of the method to treat unsteady cases was demonstrated. A situation of delayed system activation was considered. Indeed, it is to deal with such unsteady cases that the present method was developed. The results may be summed up in three steps. Firstly, ice builds up. Secondly, a system is activated (modelled by an internal heating). The method predicted the propagation of the melting front inside the ice block. Moreover, the associated temperature gradient was clearly visible. Thirdly, once the ice block at the leading edge was melted, water flowed downstream. Once out of the heated area, the water began to freeze, creating two ridges of runback ice, one on the suction side, another on the pressure side. It is worthwhile to note that Messinger's method is unable to account for such phenomena as it is based on steady state assumptions. The lack of experimental data and its high scatter does not at present time allow for a full validation of the model. However, ongoing research projects aim to fill the experimental data gap.

With these results at hand, the next step is to couple this methodology with a solver for an electrothermal ice protection system [START_REF] Chauvin | A robust coupling algorithm applied to thermal ice protection system unsteady modeling[END_REF]. Moreover, combining them with recent work on ice shedding mechanisms [START_REF] Bennani | Numerical simulation and modeling of ice shedding: Process initiation[END_REF][START_REF] Bennani | A mixed adhesion-brittle fracture model and its application to the numerical study of ice shedding mechanisms[END_REF] would enable the simulation of a full de-icing cycle. Therefore, although the x component is not zero strictly speaking, not taking it into account is a reasonable approximation. This allows to reduce the phase change problem to one dimension and to construct a model that retains the relevant physics while being much cheaper computationally.

Bibliography

  Aerodynamic tangential pressure gradient [N.m -3 ] Γ k Mass transfer term for the k-th layer [kg.m -2 .s -1 ] Tf Mean temperature of the running liquid film [K] vx Mean velocity of the running liquid film [m.s -1 ] λ k Thermal conductivity of the material in layer k [W.m -1 .K -1 ] µ w Viscosity of water [P a.s] Φ f Energy transfer term for the liquid running film [J.m -2 .s -1 ] ρ k Density of the material in layer k [kg.m -3 ] τ Aerodynamic shear stress [N.m -2 ] a k Position of the lower boundary of the k-th layer in the z direction [m] b k Position of the upper boundary of the k-th layer in the z direction [m] c k Specific heat of the material in layer k [J.kg -1 .K -1 ] g x Gravity component along the curvilinear direction [m.s -2 ]

Figure 1 :

 1 Figure 1: State of the art icing code architecture. The work presented in this paper focuses on the ice accretion solver.

Figure 2 :

 2 Figure 2: Illustration of a generic icing situation

Figure 3 :

 3 Figure 3: Illustration of the three layers: h f , h i and hs are the heights of the running film, ice and static (melted) film. The a k 's and b k 's represent the positions of the boundaries between each layer in the z direction (normal to the surface)

  (2) and (3) with h f = 0, ∂ρwh f ∂t = 0 and ∂ρwcwh f Tf ∂t = 0. The source terms for (2) and (3) are given by: Γ f e = ṁimpṁev ( Tfe ) (9) Φfe = ṁimp c w T tot imp -T f eṁev ( Tfe )L v (T f e ) + Φconv (T f e ) + Φwall (T f e ) (10) 2 -Running wet In this mode, only the liquid running film layer is present. The unknowns are the running liquid film height h f and the mean film temperature Tf . The source terms for (2) and (

Figure 5 :← l + 1 [

 51 Figure 5: Illustration of the trial process during implicit time marching. A mode is assumed and the resulting solution is checked for consistency

  (a) Comparison between the analytical solution and numerical solutions obtained with increasing number of degrees of freedom (fixed time step ∆t = 0.01s) Comparison between the analytical solution and numerical solutions obtained with different time steps (fixed DoF=8)

Figure 7 :

 7 Figure 7: Melting front position as a function of time

  (a) Temperature profile at t = 0.1s (b) Temperature profile at t = 5s

Figure 8 :

 8 Figure 8: Temperature field: comparison between analytical and numerical solutions obtained with increasing number of degrees of freedom

( a )

 a Rime ice shape (accretion time: 690s) (b) Glaze ice shape (accretion time: 360s

Figure 9 :

 9 Figure 9: comparison between MiLeS2D, MESSINGER2D and experimental data for rime and glaze cases

  Figure10: Droplet impact rate, running film and ice build up modes at t = 0.25s. The plotted thicknesses dimensionless. The maximum thicknesses are h rf,max = 1.9 µm for the running film and h i,max = 6.9 µm for the ice.

Figure 11 :

 11 Figure 11: Unsteady runback process during first instants of the simulation

  (a) Two dimensional view of the ice shape and running film.

  Thickness of ice layer and running film layer.

Figure 12 :

 12 Figure 12: Ice shape at t = 20s

Figure 13 :

 13 Figure 13: Temperature field in the ice block near the leading edge. A temperature gradient is clearly visible as the melting front propagates. The dashed line represents the melting front. The film height is scaled x50

Figure 14 :

 14 Figure 14: Temperature field and layer thicknesses during the melting process

Figure 15 :

 15 Figure15: A small region of ice is melted. The resulting liquid water runs downstream.

Figure 17 :

 17 Figure 17: Runback ice shapes: effect of refinement

Figure A. 19 :

 19 Figure A.19: Comparison between x and y components of the temperature gradient at t = 30s. Arrows represent -∇T . The dashed represents the position of the melting front.

Table 1 :

 1 Physical properties of ice used for solving the Stefan problem

	2060	4185	330000	2.1	0.6	273.15

Table 2 :

 2 Aerodynamic and atmospheric conditions for the rime and glaze icing cases. The droplet size distribution is assumed to be monodisperse

	µm)

Table 3 :

 3 Aerodynamic and atmospheric conditions for the delayed activation case. The droplet size distribution is assumed to be monodisperse7.1. Step 1: ice build up

	µm)

The sequence is then initialized using the values U n k :

On both parts (upper and lower) the iterative process proceeds from downstream. At a given cell k,

The second difficulty is to manage the appearance or disappearance of a layer and the switching between modes. In a given mode md, solving equation (45) may yield a solution U (l+1) k which is incompatible with this mode. For example, assume the cell is in a rime mode at t n and that solving equation (45) yields T (l+1) i (z = 0) > T m . This solution is clearly incompatible with the rime mode. Indeed, it means that the ice should be melting and that a liquid layer should appear. Hence the mode of the cell at iteration (l + 1) can't be rime. The mode at iteration (l + 1) should be changed and the corresponding equation (45) should be solved (and so on until one finds a compatible mode). In a more general way, the trial process of different modes is illustrated in Figure 5.

Airfoil

Figure 16: Predicted runback ice shapes after 50s (obtained with a mesh size of ∆s = 0.0001m)

Mesh convergence analysis

All the previous results obtained with a mesh size of ∆s = 0.0001m. This choice was the result of a mesh refinement analysis. Figure 17 shows a clear effect of mesh refinement on the runback shapes. Indeed, for a coarse mesh, the runback ice shapes appear to be more smoothed out. Mesh refinement has the effect of sharpening the ice shapes.

Appendix A. Neglecting the tangential heat fluxes in the ice and static film layers: approximation assessment

In the model presented previously,heat conduction is considered in the normal direction only. So as to assess this approximation, consider the case illustrated in Figure A.18 where an internal heating is applied to a specific region of an ice block. This case is simulated with an in-house code that solves the two dimensional heat transfer problem with phase change using an enthalpy method [START_REF] Bennani | Two dimensional modelling of electrothermal ice protection systems[END_REF]. This way, the temperature gradient can be computed.


Anibal Ollero
email: aollero@us.es

Guillermo Heredia

Antonio Franchi
email: antonio.franchi@laas.fr

Gianluca Antonelli
email: gianluca.antonelli@unicas.it

Konstantin Kondak
email: konstantin.kondak@dlr.de

Alberto Sanfeliu
email: sanfeliu@iri.upc.edu

Antidio Viguria
email: aviguria@catec.aero

Ramiro Martinez-De Dios

Francesco Pierri
email: francesco.pierri@unibas.it

Juan Cortés
email: juan.cortes@laas.fr

J Ramiro Martinez-De Dios

Angel Santamaria-Navarro

Miguel Angel Trujillo
email: matru-jillo@catec.aero

Ribin Balachandran
email: ribin.balachandran@dlr.de

Juan Andrade-Cetto

Angel Rodriguez

The AEROARMS Project: Aerial Robots with Advanced Manipulation Capabilities for Inspection and Maintenance

Keywords: aerial robotic manipulation, aerial robotics, inspection and maintenance

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

A ERIAL manipulation deals with aerial robots equipped with arms (see Fig. 1) or other end-effectors performing tasks, such as assembly or contact inspection, in locations that are inaccessible, very dangerous or costly to be accessed from ground. It involves very challenging problems, mainly due to aerodynamic influence and the fast dynamics when compared to other floating-base manipulation systems.

Most of the works in aerial robotic manipulation have been performed indoors, with very few exceptions such as grasping of objects from the air with a helicopter [START_REF] Pounds | Grasping from the air: Hovering capture and load stability[END_REF]. These indoor works included picking and basic manipulation with quadrotors and also cooperative assembly with several aerial robots [START_REF] Michael | Cooperative manipulation and transportation with aerial robots[END_REF]. In [START_REF] Fumagalli | Developing an aerial manipulator prototype: Physical interaction with the environment[END_REF], a quadrotor with a small arm maintaining contact with a surface was reported. Indoor redundant manipulators are studied in [START_REF] Korpela | Flight stability in aerial redundant manipulators[END_REF]. In [START_REF] Hunt | 3d printing with flying robots[END_REF], an indoor flying robot capable of depositing polyurethane foam is presented. Later on, new indoor aerial manipulators were presented and demonstrated for applications such as drawer opening [START_REF] Kim | Operating an unknown drawer using an aerial manipulator[END_REF]. Compliance was studied in [START_REF] Bartelds | Compliant aerial manipulators: Toward a new generation of aerial robotic workers[END_REF]. The transportation of bars with two quadrotors was presented in [START_REF] Nguyen | Aerial tool operation system using quadrotors as rotating thrust generators[END_REF]. More recently new aerial grasping [START_REF] Seo | Aerial grasping of cylindrical object using visual servoing based on stochastic model predictive control[END_REF] and, in general, aerial systems interacting with the environment have been presented [START_REF] Mcarthur | Design of the i-boomcopter uav for environmental interaction[END_REF]. All these systems have been demonstrated indoor. In the ARCAS project (http://www.arcasproject.eu/), indoor and outdoor aerial manipulator prototypes with both helicopters (equipped with the KUKA-DLR LWR arm) and multi-rotor systems (equipped with light arms) were demonstrated, including fully autonomous visual servoing. In AEROWORKS (http://www.aeroworks2020.eu/) collaborative exploration, aerial manipulation, and also estimation, localization and control of UAVs are included. Aerial Manipulation has been also demonstrated indoors [START_REF] Wopereis | Application of substantial and sustained force to vertical surfaces using a quadrotor[END_REF].

On the other hand, the AEROARMS project (https://aeroarms-project.eu/) is more focused on aerial manipulation for the outdoor application to Inspection and Maintenance (I&M) and it includes the development of new aerial robotic manipulators. Three use cases in industrial environments, particularly in Oil & Gas plants (see Fig. 2) were selected: direct contact measurements while flying, robotic crawler deployment, and installation of sensors. Performing experiments in these plants require previous extensive outdoor testing with reliable aerial manipulators able to deal with the harder outdoor constraints in industrial settings including the agility to react to wind perturbations, compliance to afford unexpected impacts, variations in lighting conditions, and the need to compensate for inaccuracies in GNSS positioning.

Solving the above inspection use-cases in realistic settings is a hard problem. We believe the design and development of such manipulation systems requires the following developments:

• Aerial robotic manipulators able to apply forces for contact inspection in any direction, ability to compensate wind perturbations, as well as dual arm manipulators for the installation of sensors in complex settings; • Control systems for the above aerial manipulators, integrating the aerial platform and the manipulator; • Bilateral teleoperation systems with haptic interfaces and appropriated trade-off between stability and performance; • Reliable outdoor autonomous perception system;

• Planning system considering the dynamic behavior in closed loop of the aerial platform and the arm. The above components are summarized in the following sections. The paper also presents a new robot for industrial applications, the lessons learned and the conclusions.

PLATFORMS

Aerial manipulators are multi-body systems with coupled (aerial platform-manipulator) dynamic behavior applying forces to objects. AEROARMS includes the modeling and control of propulsion systems [START_REF] Franchi | Adaptive closed-loop speed control of bldc motors with applications to multi-rotor aerial vehicles[END_REF] and aerodynamic effects when flying very close to the environment [START_REF] Sanchez-Cuevas | Characterization of the aerodynamic ground effect and its influence in multirotor control[END_REF], but this paper concentrates on the main design aspects of platforms.

Most aerial manipulators [START_REF] Fumagalli | Developing an aerial manipulator prototype: Physical interaction with the environment[END_REF]- [START_REF] Kim | Operating an unknown drawer using an aerial manipulator[END_REF] use standard multi-rotor platforms with all propellers oriented in the same direction. They are thus underactuated with only four DoF, and they need to tilt to move or exert forces in the horizontal plane. Non-conventional multi-directional thrust platforms with tilted rotors [START_REF] Ryll | 6D physical interaction with a fully actuated aerial robot[END_REF] have been designed with the capability of directing the total thrust in more than one direction in body frame (see Fig. 1(a)), thus capable of full 6D wrench exertion without tilting. This configuration has been used in the Tilt-Hex hexarotor platform (Fig. 1(a)) and in the AEROX octorotor industrial platform (Fig. 1(b)). The characteristics of these platforms are shown in Table 1.

Aerial manipulators use mechatronic devices as bars, grippers or multi-link arms to reach the operation point, grasp and manipulate objects and exert forces on the environment. The design of these devices is important since their movement induce the displacement of the Center off Mass (CoM) and variation of moments of inertia, generating disturbance torques and affecting the dynamics of the aerial vehicle [START_REF] Korpela | Flight stability in aerial redundant manipulators[END_REF]. To overcome these effects, the AEROX new arm configuration has been designed. The arms weight is compensated with the batteries, thus the movement of the arm does not modify the position of the CoM of the aerial manipulator. Additionally, any contact forces from the interaction with the environment are transmitted directly to the CoM without generating perturbation torques.

More advanced manipulation than simple grasping or exerting forces over a surface requires the use of aerial manipulators with dual arms. Very few dual arm aerial manipulators have been developed, as for example [START_REF] Korpela | Towards valve turning using a dual-arm aerial manipulator[END_REF], which use two 3-DoF arms for valve turning. In AEROARMS, two new anthropomorphic dual-arms aerial manipulators with 4 and 5 joints (see Table 1) have been developed, since this configuration maximizes dual arm manipulability. The first one is an octorotor with arms of 0.9 Kg, and the maximum load that can be lifted by the two arms with stiff joints is 0.7 kg (Fig. 1(c)). Taking into account that collisions and impacts can severely affect stability of the aerial manipulators, and compliance has been found to be very useful [START_REF] Bartelds | Compliant aerial manipulators: Toward a new generation of aerial robotic workers[END_REF], a dual arm aerial manipulator with compliant joints [START_REF] Suarez | Physical-virtual impedance control in ultralightweight and compliant dual-arm aerial manipulators[END_REF] has been developed and tested outdoor (Fig. 1(d)). It is a hexarotor and the arms have joints with elastic elements that can absorb contact and collision forces generated by physical interaction, providing a compliant behavior that increases robustness and stability and also enables measuring interaction forces for force control. The total weight for both arms is 1.3 kg and the maximum load that each arm can lift is 0.2 kg.

CONTROL SYSTEMS

One of the key advances of AEROARMS has been to deeply investigate the potentiality and limitations of control architectures for aerial manipulation. More in detail, as far as the interaction control is concerned, the solution based on a multi-directional thruster platform and a 6D flying end effector [START_REF] Ryll | 6D physical interaction with a fully actuated aerial robot[END_REF] has represented a successful integration of known robotic algorithms, dealing with both motion and interaction control, in a platform with a minimal sensor suite (pose sensor plus IMU) not even needing a force sensor. Stable interaction was obtained also in previous approaches, e.g. in [START_REF] Fumagalli | Developing an aerial manipulator prototype: Physical interaction with the environment[END_REF] where an aerial manipulator interacts with a vertical surface, but in AEROARMS there are experiments involving complex interaction tasks such as peg-in-hole with tilted holes and sliding on tilted surfaces.

Two main motion control approaches have been considered. The first one is based on modeling the aerial manipulator as an unique system and designing the control scheme on the complete kinematic and dynamic models, see for instance [START_REF] Heredia | Control of a multirotor outdoor aerial manipulator[END_REF] where an integral backstepping controller has been implemented. The second approach considers the aerial platform and the robotic arm as two separate entities and, thus, two autonomous controllers are designed, one for the aerial vehicle, where the effects of the arm on dynamics are handled as disturbances, and one for the arm, where, vice versa, the disturbance is due to the vehicle. Such paradigm has been experimentally exploited, e.g., in [START_REF] Ruggiero | A multilayer control for multirotor uavs equipped with a servo robot arm[END_REF], where the aerial vehicle controller includes an estimator of external forces and moments to compensate neglected dynamics and the arm dynamics.

Vision-based techniques have also been exploited. In [START_REF] Santamaria-Navarro | Uncalibrated visual servo for unmanned aerial manipulation[END_REF] the feedback output from a camera attached to the end effector is adopted in a hierarchical control law. The camera images are exploited to drive the arm end-effector to a desired position and orientation.

One of the AEROARMS objectives is the inspection through contact while flying. This requires to exert forces on a surface in order to maintain the sensor (i.e. ultrasonic) installed in the manipulator end-effector in contact with the surface. Since interaction forces and moments can cause large deviation in motion control and/or destabilize the system, some active and/or passive compliance should be added to the vehicle and the manipulator end-effector. Passive compliance has been adopted in the platform in Fig. 1(d), through the elastic actuators of the compliant arms. This allows to cope efficiently with interaction forces, providing also measurements of these forces that are used to close the feedback loop [START_REF] Suarez | Physical-virtual impedance control in ultralightweight and compliant dual-arm aerial manipulators[END_REF]. Adoption of active compliance approaches requires measurements of the interaction wrenches, that can be achieved, e.g., by using a wrist mounted force/torque sensor on the arm, as in [START_REF] Cataldi | Impedance control of an aerial-manipulator: Preliminary results[END_REF], where the contact wrenches are fed to an admittance filter. Such a control scheme has been tested on an underactuated quadrotor equipped with a 6-DoFs manipulator. An alternative solution to the force/torque sensor is the use of wrench estimators as proposed in [START_REF] Tomic | External wrench estimation, collision detection, and reflex reaction for flying robots[END_REF] and in [START_REF] Ryll | 6D physical interaction with a fully actuated aerial robot[END_REF]. Adoption of under-actuated platforms implies that the lateral forces in body frame, which cannot be provided by the aerial platform itself, have to be generated through the dynamical/inertial coupling between the arm and the aerial robot, which is usually very hard to be handled in an accurate way. To overcome this drawback the admittance control paradigm has been exploited by considering the platform in Fig. 1(a), characterized by noncollinear fixedly-tilted propellers which make possible to control independently the translational and angular acceleration when unconstrained, or any of the six components of the exerted wrench when in contact [START_REF] Ryll | 6D physical interaction with a fully actuated aerial robot[END_REF]. As case study has been considered the contact inspection of a pipe and a peg-in-hole emulating a sensor installation. Figure 3 shows the relevant variables for one of the experiments. Videos can be found at http://homepages.laas.fr/afranchi/robotics/?q=node/414 .

Regarding the vehicle-arms coordination, a behavioral control scheme has been proposed for handling both single and dual arm aerial manipulators. The approach in [START_REF] Muscio | Coordinated control of aerial robotic manipulators: theory and experiments[END_REF], is based on an extension of the set-based null-space-based approach [START_REF] Moe | Set-based tasks within the singularity-robust multiple task-priority inverse kinematics framework: General formulation, stability analysis and experimental results[END_REF] to the case of aerial manipulators, which involves very different challenging issues with respect to other robotic systems. From one hand, the results show that this generalization is successful, on the other hand, specific behaviors and design procedures have been tailored due to the different actuating characteristics.

TELEOPERATION

In order to evaluate the feasibility of aerial telemanipulation for industrial I&M tasks, the AEROARMS teleoperation system should ensure: 1) a stable and transparent teleoperation system under the channel characteristics of the wireless communication link, 2) stability of the coupled controller for the manipulator and aerial base.

Stable teleoperation with communication time delay

For the time-delayed teleoperation system, the trade-offs between stability and performance requirements are tackled by means of a novel four-channel architecture in which force and position signals are exchanged between both the haptic interface (master) and the aerial robot (slave). In each of the channels, a passivity-based controller is exploited to guarantee stability under nominal communication time delays, packet loss or jitter, while transparency is maintained. To this end, a widely used passivity tool, namely, time domain passivity approach (TDPA), has been utilized. The overview of the implemented bilateral controller for the AEROARMS project is shown in Fig. 4.

Stable coupled control of the manipulator-aerial base system

Stability of the coupled controller for an helicopter with a manipulator is proven in [START_REF] Kim | Passive compliance control of aerial manipulators[END_REF] for the autonomous scenario.

Here, the manipulator wrench forces are computed in the fuselage frame of the helicopter whereby passivity is ensured. As an extension to the teleoperation case, the passivity check with TDPA to remove the destabilizing effects of communication time delay is applied for the coupled controller. The onground Hardware-In-the-Loop (HIL) simulator [START_REF] Artigas | The oos-sim: An on-ground simulation facility for on-orbit servicing robotic operations[END_REF] is adapted to reproduce the dynamics and control of the aerial system to repeatedly test the bilateral control under defined conditions. Fig. 5 shows the results of the teleoperation peg-in-hole experiments performed on the HIL with the helicopter control and the dynamic simulation. Although the designed controller produced a stable and high performing system (in terms of pure teleoperation), as it can be seen, the base of the manipulator moves as a result of the manipulator motion and external interaction. It was learned that for highly intricate tasks like precise end-effector positioning and accurate force exertion on the environment, reactive dynamics of the aerial base makes the task completion highly challenging to the operator (note especially the time between 27-30 seconds in the plots). In order to aid the operator, benefits of task dependent autonomy can be made use of. So, using virtual fixtures [START_REF] Abbott | Haptic virtual fixtures for robot-assisted manipulation[END_REF] and visionbased shared control [START_REF] Abi-Farraj | A visual-based shared control architecture for remote telemanipulation[END_REF], rather than pure teleoperation are planned during the scope of AEROARMS.

PERCEPTION

Existing techniques for perception in aerial manipulation are, in most cases, methods adapted from aerial vehicles without manipulation skills to provide positional accuracy and place recognition [START_REF] Karrer | Towards globally consistent visual-inertial SLAM for collaborative aerial navigation[END_REF]. The AEROARMS I&M usecases require however, perception modules that go beyond the state of the art to: (1) accurately localize the vehicle, both during the navigation and manipulation phases; (2) localize In the AEROARMS I&M use cases, two operational stages can be distinguished: navigation in the environment and manipulation in close distance to obstacles. Real time accurate and robust 6 DoF pose estimates of the robot are computed on-board by a multi-sensor filter based on Monte Carlo localization (MCL) that integrates the measurements from the 3D LIDAR, the stereo camera and the UWB nodes using a precomputed multi-sensor map.

In flights prior to I&M missions, a multi-sensor mapping scheme builds separate maps for each sensor and treats them as different layers of the multi-sensor map, all layers sharing the same temporal and spatial frames. We preferred this scheme versus a monolithic approach because it provides higher modularity, flexibility, efficiency and gives capacity to guide the development leveraging the contribution of each sensor.

Figure 6(c) shows the multi-sensor map obtained in the experiment in Fig. 6(b), which was also the testing scenario in Fig. 10. The robot localization computed by the multi-sensor MCL in that experiment had a mean error w.r.t. RTK GPS of 9.9 cm, which was 47% lower than the robot localization estimation obtained using LOAM [START_REF] Zhang | LOAM: Lidar odometry and mapping in realtime[END_REF].

Robot pose estimation during manipulation requires higher levels of accuracy. To that end, the PL-SLAM algorithm was developed to leverage the state of the art in SLAM, ORB-SLAM [START_REF] Mur-Artal | ORB-SLAM: A versatile and accurate monocular SLAM system[END_REF], with the simultaneous estimation of point and line features [START_REF] Pumarola | PL-SLAM: Real-time monocular visual slam with points and lines[END_REF], with precise localization at a rate of 3 Hz, resulting in better accuracy than direct methods. Parameterizing lines by their endpoints provides robustness against occlusions and in poorly textured environments. To improve the localization precision a new method using camera images and Deep Learning was used, reaching a precision of 4,5cm. Fig. 7(a) shows the deep learning architecture. Faster odometry estimates are obtained by fusing IMU and optical flow data. In AEROARMS, sensor synchronization is addressed with time stamps enabling timing errors of some milliseconds, sufficient to obtain the required GNSS-free pose estimation errors.

Once the pipe is localized, the industrial use-cases require to identify and track specific I&M characteristics and artifacts on the pipe, such as welding marks or corrosion points. The specific characteristics to be tracked can be chosen on-line by an operator, thus off-line learning methods are not feasible due to the low size of the training data set. We developed an on-line semi-supervised boosting method [START_REF] Villamizar | Boosted random ferns for object detection[END_REF] that adapts to real-time changes in illumination, cast shadows or partial occlusions. In this case, the operator selects on-line the type of defect, and the perception system adapts the pattern shape, color and texture to detect and track the defect. This new method can work in two modes: in the first mode the system is as "on-line learning" and if the defect characteristics change significantly, the system asks the operator to incorporate these new characteristics in the model; and in the second mode, the system works only identifying or tracking the defect, allowing some adaptation, but big changes are not allowed. An example of weld detection and tracking using this method is shown in Fig. 6(d).

The deployment of a robot crawler is also a use case as mentioned above. Pick-up and release of the crawler by the aerial robot is performed in two steps. In the approach maneuver the crawler pose is roughly estimated with an appearance-based Deep Learning method whose architecture is shown in Fig. 7(b) and, for the actual pick up or release, a positioning system based on deformable markers was built [START_REF] Amor | On-board real-time pose estimation for uavs using deformable visual contour registration[END_REF], guaranteeing very high accuracy, 2cm precision as shown in Fig. 6(e).

MOTION PLANNING AND NAVIGATION

Motion planning is required to: 1) move the aerial manipulator from a take-off position to the area to be inspected/maintained; 2) perform the inspection/maintenance task. The former may require navigation in environments cluttered with pipes, structures, etc. (see Fig. 2). The latter requires treating dynamic constraints in a subspace defined by the task. Both methods are presented below.

Navigation in cluttered environments

A motion planning approach for aerial manipulators exploiting the concept of differential flatness has been developed.

C y C q i+1 q i y i y i+1 q(t) y d (t) Γ(q, v, x d y )
C obs Fig. 8. Schematic representation of the control-aware planning paradigm. To each trajectory in the task space corresponds a motion of the robot is generated by the control method.

Differential flatness is very useful because it allows to plan trajectories directly in a lower dimensional space without the need to consider explicitly the fulfillment of the complex nonlinear and underactuated dynamics of the system. For instance, it has been exploited for dynamic grasping using aerial manipulators with a single-link arm [START_REF] Thomas | Toward autonomous avian-inspired grasping for micro aerial vehicles[END_REF]. In the context of AEROARMS, in which an aerial manipulator has to interact by means of the end effector, it results very useful if the end-effector position and orientation belong to a flat output of the system. This is unfortunately not the case for any aerial manipulator. Our main contribution in this field has been to prove that there exists a class of aerial manipulator designs, that we named protocentric, for which this property is ensured [START_REF] Yüksel | Differential flatness and control of protocentric aerial manipulators with any number of arms and mixed rigid-/elastic-joints[END_REF].

Motion planning for task-constrained aerial manipulation

An approach that reinforces the connection between motion planning and control in the context of aerial manipulation has been developed [START_REF] Tognon | Control-aware motion planning for task-constrained aerial manipulation[END_REF]. The underlying idea consists in using the controller as a local method to connect neighboring states within a (global) motion planning algorithm. The use of the controller inside the planner guarantees the feasibility of the trajectory for the real system, satisfying task constraints in addition to geometric, kinematic and dynamic constraints. Furthermore, it also allows to better predict the behavior of the closed loop system. Another advantage is that, using a control method that directly treats the redundancy of the system, allows the planner to search for a solution directly in the reduced and more relevant task space. Planning directly at the task level enables a more straight formulation of task-constrained motion planning problems, and generally reduces the dimensionality of the search space. Finally, properly defining the task, one can choose a good trade-off between the dimension of the search space and the delegation of some DoF to the controller.

The proposed control-aware motion planner applies a sampling-based algorithm to search for a solution directly in the task state space X y = C y × V y × A y , which in general represents all possible poses (C y ), velocities (V y ) and accelerations (A y ) of the end-effector while constrained to a given task. Figure 8 shows a scheme of the approach.

The planner is able to find trajectories in constrained environments where a purely reactive gradient-based approach fails due to local minima traps. Figure 9 shows one of the scenarios used for the evaluation in simulation reproducing a plausible I&M application by contact which requires the endeffector to be always in contact and perpendicular to the pipe. The planner requires about 20 seconds to obtain a solution trajectory in this context.

AEROARMS AERIAL ROBOT FOR INDUSTRIAL CONTACT INSPECTION

Various control, perception and planning algorithms and functionalities have been described in the previous sections. After a research period where different approaches and functionalities were developed and tested in controlled environment, all this know-how has been applied to develop an aerial robot for industrial outdoor contact inspection (see Fig. 10), which is the first use case of AEROARMS.

The main characteristics of this aerial platform, called AEROX, are shown in Table I. It has a 8 rotors configuration with a small angular gap at the front and back of the robot to allow the robotic arm pass between the propellers. Their motors are tilted alternatively 30 degrees around their arms, allowing the control of the lateral accelerations to avoid the tilting of the platform. Moreover, this aerial robot has a special design integrating a new type of aerial robotic manipulator, able to rotate around the vehicle CoM, being capable of robust contact operations that can be used for precise industrial inspections. This novel robotic manipulator provides the following functionalities: orientate the end effector in one degree of freedom (up and down) with 180 working space, reject efficiently the external perturbations (wind, force and torques generated from the contact, etc.), maintain a specific force in contact with the surface and calculate the position of the end effector with respect to the aerial robot.

In order to validate this configuration, many experiments have been conducted. The video at https://youtu.be/vc6E 3iS690 shows a summary of the experiments including laboratory tests in an indoor testbed with a motion tracking localization system used as ground truth, outdoor experiments with a safety tether and outdoor free flight experiments.

Fig. 11 shows how multiple perturbations were injected using a rope in order to test its rejection and stability. Each perturbation is easily identified looking to the first graph (attitude) due to the rolling produced each time the rope is Fig. 9. Illustration of a planned trajectory for pipe inspection. The planned trajectory allows the robot to maintain the contact between the sensor (installed at the end of the arm) and the inspected pipe while avoiding collisions with the other pipe on top. The path followed by the sensor is represented by the red line. A teleoperation device, specifically designed for contact inspection operations, has been developed. It is composed of the following elements:

1) Video camera allowing remote view of the interested area; 2) Wheeled system to move the end effector when in contact with the pipe; 3) Inspection sensor: ultrasound or eddy current sensor to perform industrial contact inspections; 4) Ground Interface with camera and sensor device allowing the operator to monitor and modify all the inspection parameters. Then, the procedure to perform inspections with this novel aerial robot is:

1) The aerial robot takes-off and gets closer to the area to be inspected. This phase could be autonomous or teleoperated and may eventually benefit from the above mentioned trajectory planning; 2) The aerial robot reaches contact with the pipe to be inspected. This phase can also be autonomous or teleoperated and benefits from the perception system; 3) The aerial robot is changed into inspection mode, where the aerial robot automatically maintains its position with respect to the point of contact. From this moment, the inspection operator has full control over the system: a) The operator, looking at the camera, can move the sensor along the pipe disregarding the individual DoFs of the aerial robot, benefit from the task-priority controller, and eventually resorting to shared control methods; b) Once at the area of interest, the operator uses the sensor control interface to perform the inspection. 4) The operator can then move the end-effector to the next area of interest, or finalize the inspection; 5) After the complete inspection is finalized, the pilot switches back to flying-mode and de-attaches the aerial robot from the pipe.

LESSONS LEARNED

Experiments in industrial settings can only be performed after testing extensively the robustness of the proposed solutions. In this testing, the main lessons learned were:

• Wind and aerodynamic perturbations generated by close surfaces should be compensated by the controller fast enough to maintain safety and achieve accuracy and relative stability. Multi-directional thrust platforms are better suited to maintain forces in any direction as required for contact inspection. Compliance in the arms is also very relevant to maintain safety. Currently available light servo-controllers impose limitations in the accuracy that can be obtained in the end-effector; • The above mentioned perturbations have significant impact in navigation and manipulation, generating deviations from trajectories planned with conventional methods. The use of a controller inside the planner is useful to predict the dynamic behavior of the closed loop system. Moreover, reactivity methods are needed for dynamic environments or incremental map building;

• For time-delayed teleoperation, task-dependent shared control algorithms might be helpful to: a) track the master position accurately by the slave manipulator in spite of the base movement, and b) mask the dynamics of the aerial base to the operator so that he perceives a fixed base teleoperation; • Sensor data fusion methods should be applied to compensate the limitation of individual sensors, such as the loss of GPS, lighting changes and electromagnetic perturbations.

CONCLUSIONS

AEROARMS has demonstrated the application of aerial robots to contact inspection of industrial pipes while flying, being able to perform ultrasonic measurements to determine the thickness of the pipes wall. The results presented in this paper demonstrated the suitability of aerial robotic manipulation for these outdoor applications. Significant advances with respect to the state of the art have been necessary from design, perception, teleoperation and motion planning and control.

Most of the proposed methods and technologies have been demonstrated outdoor, which is a relevant departure from other published works on aerial robotic manipulation. The technologies of multi-directional thrusters and compliance in the end-effectors have been used in an industrial prototype for the ultrasonic measurements of pipes wall thickness. This work obtained the Overall Innovation Radar Prize 2017 of the European Commission in the ICT event (see http://www.euronews.com/2017/12/04/new-dronetechnology-wins-innovation-radar-prize-2017).

Implementation in factories of the flying-based contact measurement and the deployment of the crawler are in preparation.

Fig. 1 .

 1 Fig. 1. AEROARMS platforms: (a) Multi-directional thrust hexarotor with a rigidly attached end-effector; (b) Multi-directional thrust octorotor with an inspection arm; (c) Aerial dual-arms manipulator with stiff joints; (d) Aerial dual-arms manipulator with compliant joints

Fig. 2 .

 2 Fig. 2. Oil and Gas refinery in Germany.

Fig. 3 .

 3 Fig. 3. Experimental results for interaction control

Fig. 4 .

 4 Fig. 4. Teleoperation architecture for the aerial scenario.

Fig. 5 .

 5 Fig. 5. Position of the master and slave end-effectors (top 3) for a peg-in-hole task experiment with 40 ms round-trip delay. Notice that in the last 2 plots, the base positions and orientations change significantly with respect to the slave motion

Fig. 6 .

 6 Fig. 6. Aerial vehicle with multi-sensor setup (a). Picture during mapping experiments (b). Multi-sensor map (c) built with LIDAR -orange dots in the figure, UWB -green dots-and stereo vision. Pipe weld tracking (d). Precise localization of a robot crawler using artificial markers (e).

Fig. 7 .

 7 Fig. 7. Deep learning for precise aerial vehicle localization (a). Deep learning architecture for aerial vehicle localization (b). Deep learning architecture for crawler localization.

Fig. 10 .

 10 Fig. 10. Outdoor contact inspection experiments with AEROARMS aerial robot

Fig. 11 .

 11 Fig. 11. Indoor autonomous contact operation telemetry during an experiment where 17 strong perturbations up to 100 N were injected. The attitude graph (top) shows the low effect of the perturbations on the robot attitude. The robot 3D position estimated error in the experiment (bottom) was in average lower than 1 cm even with these perturbations.

Preprint version final version at http://ieeexplore.ieee.org/ 2 Accepted for IEEE Robotics and Automation Magazine 2018

Preprint version final version at http://ieeexplore.ieee.org/ 3 Accepted for IEEE Robotics and Automation Magazine 2018

Preprint version final version at http://ieeexplore.ieee.org/ 4 Accepted for IEEE Robotics and Automation Magazine 2018

Preprint version final version at http://ieeexplore.ieee.org/ 5 Accepted for IEEE Robotics and Automation Magazine 2018

Preprint version final version at http://ieeexplore.ieee.org/

Preprint version final version at http://ieeexplore.ieee.org/ 7 Accepted for IEEE Robotics and Automation Magazine 2018

Preprint version final version at http://ieeexplore.ieee.org/ 8 Accepted for IEEE Robotics and Automation Magazine 2018

Preprint version final version at http://ieeexplore.ieee.org/ 9 Accepted for IEEE Robotics and Automation Magazine 2018

ACKNOWLEDGMENTS

The AEROARMS Consortium consists of 10 European partners (https://aeroarms-project.eu/) led by the Univ. Seville. The funding of the European Commission under the H2020 contract 644271 is acknowledged.