Proposition of a modal filtering method to enhance heat source computation within heterogeneous thermomechanical problems

Keywords: Thermal field, Infrared thermography, Spectral projection, Modal decomposition, Heat sources estimation

This paper presents a new approach to evaluate heat sources from thermal field measurements. A modal projection based on dynamics (Discrete Modal Decomposition) is used to estimate the spatial term of a heat diffusion problem. A numerical example is presented in order to validate this approach and compare it to a more classical spectral decomposition (based on thermal considerations). Results show that the proposed projection basis not only provides closer assessment of the heat sources but is also more stable to noise and side effects. Finally, a basis enrichement method is presented and tested. It shows that a priori knowledge of the heat sources shape though approximated (e.g. from strain measurement) enhances the assessement of calorific effects accompanying material deformation.

Introduction

The study of calorific effects associated with material deformation has been widely investigated since the early works of Lord Kelvin. In the late 90s, experimenters have turned to the use of infrared (IR) thermography in order to estimate the temperature variations at the sample surface during mechanical loading. The observed temperature variations depend on both the material and the type of loading and can, under some assumptions, provide valuable data on the material behaviour. For example, positive temperature variations in one area of a loaded sample may reveal that irreversibilities take place at this location. However, the intensity of the temperature variations does not only depend on the internal mechanisms of the material (reversible or not). They also depends on various characteristics of the thermal problem (geometry of the studied sample, thermophysical properties, boundary and initial conditions). Within this context, another approach, based on heat sources analysis, was proposed and has proved its worth on various thermomechanical problems [START_REF] Chrysochoos | Plastic and dissipated work and stored energy[END_REF][START_REF] Chrysochoos | An infrared image processing to analyse the calorific effects accompanying strain localisation[END_REF][START_REF] Louche | Thermal and dissipative effects accompanying Lü d e r s b a n d p r o p a g a t i o n , M a t e r i a l s S c i e n c e a n d E n g i n e e r i n g[END_REF][START_REF] Louche | Thermal observations associated with the Portevin-Le Chatelier effect in an Al-Mga l l o y ,M a t e r i a l sS c[END_REF].

Hence, heat sources has become relevant quantities to study material behaviour. They provide complementary informations on the instantaneous energy, locally associated with the material behaviour. Either the sign and the magnitude of this volumic power (denoted w 0 ch )c a nb el i n k e dt ot h e material underlying physics and can also be used to achieve energy balance [START_REF] Chrysochoos | Fields of stored energy associated with localized necking of steel[END_REF][START_REF] Dumoulin | Heat sources, energy storage and dissipation in high-strength steels: Experiments and modelling[END_REF].

Classically, the heat sources computation relies on performing both spatial and time derivation of the measured temperature fields T (x, y, t). Under some classical assumptions [START_REF] Chrysochoos | An infrared image processing to analyse the calorific effects accompanying strain localisation[END_REF], the 2D-heat diffusion equation is given by

ρC ∂θ ∂t -k ✓ ∂ 2 θ ∂x 2 + ∂ 2 θ ∂y 2 ◆ = w 0 ch with θ(x, y, t)=T -T 0 , (1)
where θ stands for the temperature variation, x and y are the spatial coordinates, ρ is the mass density, C the thermal capacity, k the thermal conductivity, T 0 (x, y)theinitialtemperaturefieldandw 0 ch the internal heat sources. Such a definition allows the assessment of w 0 ch from the time-derivative and the second space-derivative of θ(x, y, t)a n dt h ea priori knowledge of the thermophysical constants ρ, C and k.

In practice, the time resolution of modern IR cameras (up to 380 Hz) often allows the computation of a proper time derivative of θ using low-pass or polynomial filters at every pixel. On the other hand, the space-derivative term cannot be estimated as straight-forwardly. The measurement noise, intrinsic to IR-thermography, is amplified by space-derivation operators and massively impairs the calculation of the Laplace terms in Eq. [START_REF] Chrysochoos | Plastic and dissipated work and stored energy[END_REF].

Many numerical approaches have been developped in the last decade in order to estimate as accurately as possible the left hand side term of Eq. ( 1). These regularisation methods can be classified into 3 main categories:

• Filtering methods, aimed at supressing the noise from the measurement using Gaussian [START_REF] Xiaoyan | Improving the peak wavelength detection accuracy of Sn-doped, H2-loaded FBG high temperature sensors by wavelet filter and Gaussian curve fitting[END_REF][START_REF] Bozzoli | Experimental validation of the filtering technique approach applied to the restoration of the heat source field[END_REF] or low pass filters based on linear, Fourier or wavelet definitions [START_REF] Chrysochoos | An infrared image processing to analyse the calorific effects accompanying strain localisation[END_REF][START_REF] Yi | Source terms identification for the diffusion equation[END_REF][START_REF] Chuli | Wavelet and error estimation of surface heat flux[END_REF].

• Inverse methods, aimed at obtaining a quasi-solution through the use of one or several optimization algorithm [START_REF] Wong | A hybrid numerical/experimental technique for determining the heat dissipated during low cycle fatigue[END_REF][START_REF] Kaipio | Statistical inverse problems : Discretization, model reduction and inverse crimes[END_REF].

• Projection methods, they are among classical means in inverse heat transfer analysis. They consist of decomposing the measurement within a spectral basis built from eigen function of the Laplace operator [START_REF] Chrysochoos | An infrared image processing to analyse the calorific effects accompanying strain localisation[END_REF][START_REF] Doudard | A probabilistic twoscale model for high cycle fatigue life predictions[END_REF]. These techniques have recently experienced various improvements through the use of wavelet decompositions [START_REF] Candau | Solving inverse heat diffusion problem with discrete wavelet transform[END_REF]o rb r a n c hm o d e sd ecomposition [START_REF] Neveu | Simulation de la conduction non linéaire en régime variable: décomposition sur les modes de branche[END_REF]. This latter approach allows the use of generalized boundary conditions which usually represents the main disadvantage of projection methods. Mixed temporal (Fourier) and spatial (Laplace) projection methods have also been used and provided promising results [START_REF] Renault | A spectral method for the estimation of a thermomechanical heat source from infrared temperature measurements[END_REF]. Nevertheless, boundary conditions and heat sources reconstruction on the domain contour remains a significant drawback of such approaches.

The present paper proposes to implement a projection method and compare the results obtained using a decomposition built from structural dynamics and another from thermomechanical considerations. The use of a basis, not related to the solution of a given problem, is nowadays among classical means in noise removal approaches [START_REF] Adragna | Analysis of shape deviations of measured geometries with a modal basis[END_REF][START_REF] Le Goic | Multi scale modal decomposition of primary form, waviness and roughness of surfaces[END_REF][START_REF] Wang | Shape features and finite element model updating from full-field strain data[END_REF]. For the sake of comparison, an u m e r i c a le x a m p l ei sb u i l t( s e c t i o n3). The heat sources are then rebuilt using (i) the proposed modal basis (section 2.2)and(ii) the classical spectral basis (section 2.3). Hence, the influence of various parameters along with the measurement noise are investigated (section 4). Finally, a basis enrichement approach is proposed and tested (section 5).

Setup of projection bases

Projection operator

The thermal field measurement, such as provided by an IR-camera, is a film made of 2-dimensional arrays. The finite resolution of the capture leads to define an integer grid

M = R δ ⇥ R δ \ [0; a] ⇥ [0; b]whic hstepisdenotedδ.
In addition, the constant frame rate results in the measured field being part of a monotonous discrete sequence:

θ(x, y)=(θ 1 (x, y); θ 2 (x, y); ...; θ k (x, y)) , (2)
where each time step (denoted using subscript k 2 N ⇤ )i sat e m p e r a t u r e variation frame.

In a space continuous framework, the use of the modal/spectral approach requires the projection of the measured fields θ k (x, y) 2 R into an eigen basis B ⇢C 2 ([0; a] ⇥ [0; b], R) independant of the steptime k.T h i sb a s i si sd e fi n e d by its continuous eigen vectors denoted (Q c i ) i2N * . Hence, assuming that B is a basis of the solution subspace S of Eq. ( 1), the measured field can be expressed as a linear combination of the eigen vectors:

θ c k (x, y)= +1 X i=1 α c k,i • Q c i (x, y), 8(x, y) 2 R 2 (3)
However, assuming a truncation N of the sum in Eq. ( 3), a residual terms R N (x, y)a r i s e si nt h es e c o n dm e m b e r . T h i sl a t t e rt e n d st o w a r dz e r oo nM if the grid step δ tends toward zero. In other words it is easy to settle that:

lemma 1 -R N (x, y, δ)-! 0, 8(x, y) 2 R 2 i) if δ-! 0 ii) and if N = |M|
where |•| stands for the cardinality of •. This results prooves that the space continuous and discrete cases merges when the truncation N equals the number of nodes in M. Hence, if B is a basis of the solution subspace S if Eq. (1)t h e nt h ep r oj e c t i o no fθ c k on M is exact. Moreover, one may notice that for a given discretisation δ, the independance of the discrete modal vectors

Q i is sufficient to ensure that lemma 2 -R N,δ (x, y)=0, 8(x, y) 2M i) if B is a basis of the set of function of M in R, ii) and if N = |M|.
A consequence of lemma 2 is that if B is a free familly (i.e. if N< |M|), the residual term is a continuously decreasing function of N . Another consequence is that B doesn't have to be a basis of S, indeed any independent set of vectors verifies lemma 2. Subsequently, for a given space discretization, if the number of modes equals the number of nodes then the residual can be zero regardless of whether B is or isn't a basis of S. This consideration leads to assume that the relation beween the projection basis and the underlying physics is not relevant if N is high enough.

Modal basis related to dynamics

Let U be a compact subspace of R 3 and q : U ! R a C 2 -continuous function. A rectangular domain U =[ 0 ;a] ⇥ [0; b] ⇥ [t k ; t k+1 ]d e s i g n a t e st h e solid structure under investigation. The associated geometry (a square plate in the proposed problem) leads to the definition of the following dynamic structural problem:

M • q + K • q =0 with q = q(x, y, t), (4)
where M and K stand for the mass and the stiffness matrices respectively. Under such formalism, q(x, y, t)s t a n d sf o rt h ed i s p l a c e m e n t sw h i c h caracterizes the modal shapes. Such a problem classically gives a frequency based solution

q(x, y, t)= +1 X i=1 Q i (x, y) • cos(ω i t), (5)
where Q i is the magnitude vector associated with the pulsation ω i . Hence, in the finite dimension framework, the eigen modes defined by (Q i ,ω i )a r e determined by solving the following linear system

✓ M -1 K - 1 ω 2 i I ◆ Q i =0, ( 6
)
where I is the identity matrix and M -1 K is of course assumed to be diagonalizable. The discrete solution is computed by using Finite Element Analysis (FEA) and provides the dynamic modal basis

B D =(Q 1 ,Q 2 ,...,Q N )
The numerical finite element problem is solved using Abaqus software. The procedure can be summarized as following:

1. As q u a r ep l a t ei sm od e l l e da sa3 Dd e f o r m a b l es h e l l .

2. The structure is meshed using 100⇥100 square quadratic shell elements (9 nodes, normal z).

3. Boundary conditions are prescribed and impose no in-plane displacement at every nodes, the scalar product of the displacement vector q(x, y) with axis x and y being then: hq(x, y, t); xi =0andhq(x, y, t); yi = 0e v e r y w h e r e .

4. The frequency solver is run with a requested number of mode of N = 500 (CPU time is 1266s).

5. Modal displacements along z are exported to matlab and interpolated over a 300 ⇥ 300 monotonic grid (The 25 first modal vectors are ploted in Fig. 1-(a)).

However, the non orthonormality of B D does not allow the use of the classical projector P Q = QQ T .I n d e e d , t h e u s e o f t h e d u a l b a s i s

Q ⇤ = Q T Q -1 Q T is required.
An infinite norm is given to the modal vectors such as kQ i k 1 =1 . T h u s ,t h es e to fm o d a lc o o r d i n a t e sα k resulting from the projection of the displacement field at step k (denoted q k (x, y)) within a non-orthonormal basis is given as follows:

α k = {α k,i 2 R} = Q T Q -1 Q T • q k (7)
where Q denotes the matrix whose columns are Q 1 ,...,Q N ,a n dα k = (α k,i ) i2[1;N ] is the set of modal coordinates. The point of the present paper is to associate the displacement fields q(x, y, t k )t ot h em e a s u r e dt e m p e r atures θ k (x, y). Indeed, a consequence of Eq. ( 6)i st h a ti fN = |M|,t h e n B D is completed and therefore defines a basis in the solution subspace of Eq. (1)o nM. Hence, even though B D is far from describing the actual underlaying physics, it may (in theory) lead to an exact solution at nodes if N = |M|.I np r a c t i c e ,t h i sa m o u n to fm o d e s( e . g .1 0 0 0 0f o ra1 0 0 ⇥ 100 mesh) is of course not relevant from a numerical standpoint, however the obtained basis exhibits several significant advantages. First of all, it naturally describes a space vector and thus allows a straight forward projection. In addition, the energy sorting of modes leads to an increasing complexity of shape (and frequency) when the truncation order N increases. Finally, it does not require the use of variable separation and therefore leads to a slower frequency increasing along N .I no t h e rw o r d s ,t h ek th mode of B D exhibits a smaller frequency than its counterpart in the spectral basis (calculated assuming variable separation). This latter point is of the utmost interest since the computation of the second order space-derivative (finite differences) is well known for significant divergence when high frequencies are involved.

Spectral basis related to the Laplace operator

Let us consider the two dimensional thermal equation in the rectangle

U = [0; a] ⇥ [0; b] ⇥ [t k ; t k+1 ]d e fi
8 > > > > > > < > > > > > > : ∆θ = ∂ 2 θ ∂x 2 + ∂ 2 θ ∂y 2 = ρC k ∂θ ∂t over U ∂θ ∂x (0; y, t)= ∂θ ∂x (a; y, t)=0 o v er∂U x ∂θ ∂y (x;0,t)= ∂θ ∂y (x; b, t)=0 o v er∂U y , (8)
where ∂U x [ ∂U y = ∂U. Hence, if θ is solution of this problem, it is an harmonic function. Solutions are found assuming separation of variables. Thus, θ(x, y, t)=f (t)u(x)v(y). Substituting this in Eq. ( 8)a n dd i v i d i n gb y θ gives:

ρC k ḟ f = ✓ u 00 u + v 00 v ◆ . (9)
However, the independence of u, v and t to their respective variables imposes those three ratios to be constant. Then, variable separation states the problem as three independant differential equations:

u 00 = λ 1 u and v 00 = λ 2 v and ḟ = γf, (10)
such as λ 1 + λ 2 = γ.T h es h a p eo ft h es o l u t i o n s ,f o ru and v,d e pe n d so n the sign of the λ i which requires the examination of the two following cases:

λ i  0 ) 8 < : u(x)=A • cos ⇣ p |λ 1 |x ⌘ + B • sin ⇣ p |λ 1 |x ⌘ v(y)=C • cos ⇣ p |λ 2 |y ⌘ + D • sin ⇣ p |λ 2 |y ⌘ (11)
or

λ i > 0 ) 8 < : u(x)= Ã • cosh ⇣ p |λ 1 |x ⌘ + B • sinh ⇣ p |λ 1 |x ⌘ v(y)= C • cosh ⇣ p |λ 2 |y ⌘ + D • sinh ⇣ p |λ 2 |y ⌘ . (12)
Therefore, the boundary conditions of Eq. ( 8)c a nb eu s e dt od e t e r m i n e the value of the 8 constants A λ ,B λ ,..., Dλ .I tc o m e s ,t h a tt h es o l u t i o n s postulated in Eq. ( 12)l e a d st ou(x)=v(y)=08(x; y) 2 R 2 .T h u s ,t h eo n l y solutions postulated in Eq. ( 11)canbesolutionoftheboundaryproblem. In such case, it is easy to show that the boundary conditions impose:

8 > > > > < > > > > : B = D =0 λ 1 = - ⇣ nπ a ⌘ 2 and λ 2 = - ⇣ mπ b ⌘ 2 γ = - ✓ n 2 π 2 a 2 + m 2 π 2 b 2 ◆ , (13)
with n and m in Z. According to this expression of γ,t h et e m p o r a l solution is given by:

f (t)=e -σn,mt with σ n,m = kπ 2 ρC ✓ n 2 a 2 + m 2 b 2 ◆ . (14)
Then, assuming a steady state case within [t k ; t k+1 ], it can be shown that f (t) can be approximated by 1. Hence

θ k (x, y)= N X n=1 N X m=1 A n B m cos ⇣ nπx a ⌘ cos ⇣ mπy b ⌘ | {z } mode vectors Qn,m . (15)
Finally, the vectors of the spectral basis

B ∆ =(Q 1,1 ,Q 2,1 ,...,Q n,m
)associated with the Laplace operator correspond to the eigenfunctions estimated from the right hand side of Eq. [START_REF] Neveu | Simulation de la conduction non linéaire en régime variable: décomposition sur les modes de branche[END_REF]. Since the vectors are normalized, B ∆ is classically a completed and orthonormal basis on M [START_REF] Boas | Mathematical Methods in the Physical Sciences[END_REF]:

kQ 1,1 k 1 = kQ 2,1 k 1 = ... = kQ n,m k 1 =1. ( 16
)
In order to sort the obtained basis in increasing order of frequencies, the modes are sorted by increasing products n⇥m.D o u b l eo c c u r e n c e sa r et r e a t e d by placing the greater value of n ahead. One of the main consequence of the sorting operation is that it must be perfomed again if N (denoting the total number of modes or truncation order) varies. This disadvantage prevents from a single, once and for all calculation of B ∆ . Another disadvantage comes from the periodicity of the Q i vectors which requires high truncation order to discribe a non-periodic shape. However, this basis exhibits a major advantage: it is a parabolic solution of Eq. ( 1)a n di st h e r e f o r ec l o s e -t ophysics.

Modes derivation

As shown in Eq. ( 1), the measured field must be derived twice in space in order to calculate the Laplace operator. However, since the projector P Q is a linear transformation, the additivity and homogeneity properties allows the writing of the Laplace operator as:

∆θ k (x, y)= N X i=1 α k,i ✓ ∂ 2 Q i (x, y) ∂x 2 + ∂ 2 Q i (x, y) ∂y 2 ◆ +∆R n,k (x, y). ( 17
)
The point is then to assess the second derivative of the mode vectors Q i whether they are stemmed from the modal or the spectral basis. In practice, these basis are evaluated using different numerical means: B D is calculated through a finite element interpolation while B ∆ presents an analytical relation Eq. ( 15)t h a ti sf o r w a r dt oi m p l e m e n tw i t h i nM a t l a b c for example.

Modal basis

The modal vectors are calculated using a Finite Element model made of 10000 quadratic 9-node elements. This leads to an obvious continuity problem since Q i 2 B D is C 2 -continuous over each element Ω ⇢ R 2 but only C 1 over ∂Ω, which imposes ∆Q i to be a piecewise linear function C 0 . Such a property results in an inaccurate interpolation of the Q i .T h e r e f o r e , a concatenation of 4 quadratic elements is used to settle higher order interpolation functions that allows a C 2 -continuous definition of the Laplace operator. In practice, the element Ωisdefinedas Ω= S 4 i=1 Ω i . Hence, Ωisa 25-nodes square element provided with 4 th order interpolation functions (C 4continuous). It results in ∆Q i to be twice continuous over

U =[0;a] ⇥ [0; b] [18].

Spectral basis

The C 1 -continuity of the triginometric functions along with the variable separation performed within the resolution of problem Eq. ( 8)l e a d st oa straight forward analytical expression of the second derivative of

Q n,m 2 B ∆ as: ∆Q n,m (x, y)=C n,m h cos ⇣ nπx a ⌘ cos ⇣ mπy b ⌘i . ( 18
)
where

C n,m = - ⇣ nπ a 2 + mπ b 2 ⌘
is a multiplicative constant.

Validation from a numerical example

Example description

The comparison of the two presented projection bases is performed by an u m e r i c a le x a m p l e . As q u a r ep l a t ee x h i b i t i n gt w od i s t i n c th e a ts o u r c e s is modeled in Abaqus c software. The finite element diffusion problem is solved numerically over a time period of t =2 0 0s.T h er e s u l t i n gt e m p e r ature fields are output to a Matlab c environement and interpolated over the 300⇥300 monotonic grid. A random white noise is added (standard deviation is denoted σ)a n dt h ep r o j e c t o rE q .( 7) is used to approximate the original temperature fields. Finally the space derivation of modes Eq. ( 17)an dan umerical time derivation (finite differences) of the obtained temperature fields allows the use of the heat diffusion equation Eq. ( 1) and an approximation of the imposed heat sources.

The chosen geometry is a 3D plane-parallel plate with dimensions (a; a; e). As shown in Fig. 2,t w oz o n e sa r ec h o s e nt obet h ea c t i v eh e a ts o u r c e so ft h e model. The first one is rectangular and host a pyramid-shaped distribution of heat (region no.1). The second one is circular and exhibits a peak-shaped

Step #1

Diffusion Interpolation

Step #2

Step #3

Step #6

Step #5

Step #4 distribution of heat (region no.2). The boundary conditions over the six sides of the model imposes a zero output heat flux such as:

8 > > > > > > < > > > > > > : ∂θ k ∂x (0; y; z)= ∂θ k ∂x (a; y; z)=0 o v er∂U x ∂θ k ∂y (x;0;z)= ∂θ k ∂y (x; a; z)=0 o v er∂U y ∂θ k ∂z (x; y;0) = ∂θ k ∂z (x; y; c)=0 o v er∂U z , (19)
where ∂U x [ ∂U y [ ∂U z = ∂U.T h ea c t u a ld i m e n s i o n sa r ea =1 0 0mm and the thickness is e =5mm. As shown in Fig. 5,t h em a g n i t u d e so ft h e heat sources in region no.1 and no.2 evolve independently with time. Finally, the computed temperature fields are output every 10s (1  k  20) on the plate surface (z =0).

The use of a 300 ⇥ 300 monotonic grid (denoted {x i ; y j } with δ = x ix i-1 = y i -y i-1 , 8(i, j) 2 [2; 300] ⇥ [2; 300]) relies on the difference between the mesh used for mode calculation: M = {u i ; v j } and the one used in solving the thermal FE problem: {X i ; Y j }. As shown in Fig. 3,t h i sl a t t e ri sl o c a l l y refined and non-monotonic, hence the interpolation of the temperature fields over {x i ; y j } assume the role of the IR camera spatial discretization (with a similar resolution). The conductivity is arbitraly set to k =2 5W.m -1 .K -1 ,t h et h e r m a l capacity to C =520J.kg -1 .K -1 and the mass density is ρ =7800kg.m -3 .

Results: Heat sources comparison

Considering a noiseless (σ =0 )m e a s u r e dfi e l dθ k (x, y), the heat sources reconstruction is performed as presented above. The truncation order is set to N =5 0 0f o rb o t hm o d a la n ds p e c t r a lp r o j e c t i o n s . T h eo b t a i n e dh e a t sources are given in Fig. 4.T h er o o to fm e a ns q u a r ee r r o r( d e n o t e dRMS) is computed as a comparison between the imposed heat sources and the computed ones (denoted w 0 ch and w0 ch respectively):

RMS k = 300 X i=1 300 X j=1 q w0 ch,k (x i ,y j ) -w 0 ch,k (x i ,y j ) 2 . ( 20
)
It can be seen that both projections provide similar heat source reconstructions. However, the RMS errors show that the projection within B D provide slightly closer heat sources. The time history at point 1 and 2 is plotted in Fig. 5.I m p o s e dh e a ts o u r c e sa r ep l o t t e dv e r s u st h er e c o n s t r u c t e do n e s (using both B D and B ∆ projection). It is noticable that the computed heat sources do not significantly vary when one projection is used or the other. Indeed, one should conclude from the presented results that, in absence of noise, modal and spectral projections provide quasi-identical heat sources at every locations and step times.

Discussion

Sensitivity to measurement noise

The measurement noise generated by IR-camera capture may lead to significantly small signal-to-noise ratio if the recorded phenomenon involves small temperature rises. Indeed, the difficulties encountered in inverse heat transfer analysis mainly rely on the measurement noise which leads to a illposed problem which becomes very unstable to the double space derivation operator. This specific point has been adressed in many published studies [START_REF] Le Niliot | Infrared thermography applied to the resolution of inverse heat conduction problems: recovery of heat line sources and boundary conditions[END_REF][START_REF] Louche | Thermal and dissipative effects accompanying Lü d e r s b a n d p r o p a g a t i o n , M a t e r i a l s S c i e n c e a n d E n g i n e e r i n g[END_REF][START_REF] Boulanger | Calorimetric analysis of dissipative and thermoelastic effects associated with the fatigue behavior of steels[END_REF], hence the stability of the presented approach to noise is therefore highly relevant. 6: Heat sources reconstruction using modal and spectral projection (red lines) for increasing noises compared with the imposed sources (black lines). Graphs depicte the evolution of the reconstructed heat sources along two lines that cross the main heat sources (region no.1 and no.2 presented in Fig. 2).

It has been previously shown that when the noise tends toward zero, both the projection bases provide a good approximation of the heat sources. On the other hand, the increasing of the noise magnitude leads to impair the calculations. It is therefore assumed that reconstruction error exhibits a strong dependence on the signal-to-noise ratio (denoted SNR). It is classically defined as:

SNR = max |θ k (x i ,y j )| ijk σ (21)
Accordingly, a random white noise was artificially added to the numerical temperature fields with 30 increasing values of the SNR. Fig. 6 shows the reconstructed heat sources profiles along the two lines that cross regions no.1 and no.2. Moreover, the RMS errors are computed and plotted in Fig. 7. Two conclusions can be drawn from Fig. 6 and Fig. 7. Firstly, as expected, the reconstruction error increases rapidly when the SNR decreases. Secondly, the modal projection presents an improvement from the spectral approach since the RMS error is smaller for every noise magnitude. Though, the gap between these two reconstructions becomes particularly visible for higly noised measurements (see Fig. 7 left hand side), the relative improvement is actually about 10% for all the studied SNR.

T h i ss p e c i fi cp o i n tc o n s t i t u e s am a j o ri m p r o v e m e n to ft h ep r o p o s e dp r o j e c t i o nm e t h o ds i n c ei ta l l o w sa

better heat source assessment when temperatures rises are small: e.g. small strain rates, high conductivity, early localization, vicinity of peaks... Another relevant point can be seen in Fig. 6-(c). In this case, the heat sources in region no.1 (pyramid-shaped) is actually very poorly estimated while the heat sources in region no.2 is far better approximated (especially with B D ). This emphasizes the shape influence on the reconstruction capability. It is well known that square-like shape leads to the use of high frequency modes and therefore impairs the Laplacian assessment. However, as proposed in [START_REF] Adragna | Analysis of shape deviations of measured geometries with a modal basis[END_REF][START_REF] Doudard | Identification of heat source fields from infrared thermography: Determination of 'self-heating' in a dual-phase steel by using a dog bone sample[END_REF], basis enrichement can be helpfull in such case (see section 5).

Sensitivity to truncation order

The results presented above are obtained for a truncation order of N = 500. However, this choise is arbitrary, and other values lead to different approximations of the imposed heat sources. In absence of noise, it can be shown that increasing N always provides a closer approximation of the heat sources. Such an asymptotic problem imposes RMS to tend toward 0 when N tends toward the number of measured points (see section 2.1). One must notice that the use of interpolation (see Fig. 3) prevents in practice from an exact solution of the problem. However, under the assumption of σ = 0, RMS(N )i sac o n t i n u o u s l yd e c r e a s i n gf u n c t i o n . T h ei l l -p o s e dp r o b l e m , settled by the addition of noise, no longer verifies this property. Therefore, the existence of an optimal truncation order (denoted N ⇤ )c a nbea s s u m e d .

Such observations make it necessary to study, not only the influence of N on the reconstruction quality, but also the influence of the signal-to-noise ratio on N ⇤ .I np r a c t i c e ,4 0v a l u e so fN were selected beween 25 (assumed here to be the smallest reasonable number of modes for a proper reconstruction) and 500 (assumed to be sufficent to obtain a very good approximation of the measured temperature field). Beside, the influence of 10 increasing values of noise on N ⇤ are also studied.

It can be seen from Fig. 8-(a) that the modal projection provides better heat source approximations regardless of the measurement noise. As expected, in absence of noise, RMS(N )c o n t i n u o u s l yd e c r e a s ea n dt h em o d a l projection slightly prevails in such a case. One may also notice that increasing the noise magnitude leads to a bigger gap between the spectral and modal projection approaches. This tends to confirm the previous observation (section 4.1): the proposed approach improves the stability to measurement noise. Fig. 8-(b), shows that the optimum truncation order N ⇤ is always higher for the modal projection. Indeed the lower frequencies of the dynamic modes allows the use of higher truncature order even though the signal-tonoise ratio is poor, and thus provide a better temperature and heat source approximation.

Sensitivity to convective exchanges

The numerical example such as defined above does not include any dependency to convective exchanges over the main surface. Experimental ob- servations classicaly lead to consider this assumption as strong and therefore requires to evaluate the sensitivity of the presented results to convective boundary conditions. Indeed, real measurement situation often lead to temperature field smoothened by convective losses over the domain borders. In such case the heat diffusion equation Eq. ( 22)be c o m e s :

ρC ✓ ∂θ ∂t + θ τ ◆ -k ✓ ∂ 2 θ ∂x 2 + ∂ 2 θ ∂y 2 ◆ = w 0 ch ( 22
)
where τ is the time constant of the convective exchanges related to the convection coefficient h z in the z direction as following [START_REF] Chrysochoos | An infrared image processing to analyse the calorific effects accompanying strain localisation[END_REF]:

h z = ρCe 2τ ( 23
)
where e is the thickness of the plate. The value of h z was set to twelve increasing values inbetween 0 and 1000 W.m -2 .K -1 and the RMS error between the left hand side of Eq. ( 22) (computed) and the right hand side of it (imposed) was monitored. The inertial and convective terms of Eq. ( 22) are evaluated from θ and its time derivative for both spectral and modal reconstruction. The laplacian term is assessed by the mean of either spectral or modal projection. In addition, one should notice that adding convective exchanges over the main surface does not affect the definition of the spectral basis such as presented in Eq. [START_REF] Neveu | Simulation de la conduction non linéaire en régime variable: décomposition sur les modes de branche[END_REF]. Indeed the θ/τ term clearly becomes a constant when accounting for the variable separation performed between Eq. ( 8)andEq.(9)(dividingbyθ). Hence, such a formulation of the problem affects the value of the σ n,m constant but leave the modes shapes unchanged. It can be seen that for all investigated cases, the modal projection results in a lower RMS error. However the gap between the quality of the reconstructions appears to be bigger when noise is added which is consistent with the observation presented in section 4.1.I ti sa l s ov i s i b l et h a tt h eR M S error decreases when the convection magnitude increases. Indeed, a higher convection leads to smoothen the temperature fields and therefore decrease the influence of the Laplacian term of Eq. ( 22). Hence, it seems consistent to obtain a better fit with smothen fields than sharp fields.

Boundary conditions and side effects

Another common application of the thermal measurement is the estimation of the heat flux governed by convective exchange phenomenon. The Fourier law is classically used to mo del such phenomenon [START_REF] Chrysochoos | An infrared image processing to analyse the calorific effects accompanying strain localisation[END_REF]. It offers an easy way of estimating the heat transfer coefficient h from the following uni-directional solution:

-kA ∂θ ∂x i ∂Ux i = h x • A (T (t) -T env ) , (24)
where A stands for the surface area of the heat being transferred and T env is the environement temperature (here chosen to equal T 0 such as T (t) -T env = θ). Hence, one may notice that assuming the space-derivability of the measured temperature field, h x can be assessed in a straight forward manner as:

h x = - k θ ∂θ ∂x x=a . ( 25
)
As depicted in Fig. 10-(a) two convective exchanges are added to the finite element model definition and lead to thermal losses through ∂U x and ∂U y . The imposed magnitudes are h x =100W.m -2 .K -1 and h y =50W.m -2 .K -1 over two cross section areas of 100 ⇥ 5 mm 2 . All other parameters remain unchanged. As pointed by [START_REF] Chrysochoos | An infrared image processing to analyse the calorific effects accompanying strain localisation[END_REF][START_REF] Neveu | Simulation de la conduction non linéaire en régime variable: décomposition sur les modes de branche[END_REF], one of the main disadvantage of the spectral decomposition is the bad quality of the approximation at the borders. Local divergences of the RMS error arises (see Fig. 10-(b)) due to the definition of the spectral basis used in Eq. ( 8). Indeed, B ∆ is built under the assump-tion of Neumann's conditions and is thus not suited to model convective exchanges. In addition, it must be emphasized that, without a priori knowledge of the heat transfer coefficients h x and h y ,i ti spo i n t l e s st or e -s e t t l et h e Laplace problem, Eq. ( 8)i no r d e rt ob u i l das pe c t r a lb a s i st h a ta c c o u n t sf o r either Dirichlet's or mixed boundary conditions. Indeed, the fact that B ∆ is close-to-physics become a major problem when this physics is only partially known.

h y =50 W.m -2 .K -1 h x =100 W.m -2 .K -1
It can be seen in Fig. 11 that using the spectral projection disables any straight identification of the heat transfer coefficients. However, the use of a linear approximation over the vicinity of the borderline allows an approximation of h x and h y . Fig. 11 shows the results obtained for a linear interpolation from x i = a/10 to x i = a (see Fig. 10-(b)).

On the other hand, the use of modal projection provides good approximations of the imposed heat transfer coefficients and apprears, once again, to be less noise-sensitive (Fig. 11). The modal basis is build from dynamics considerations and is therefore not as close-to-physic as the sprectral one. However, the imposed boundary conditions (free translation along z and free rotation around x and y)r e s u l t si nt h ea b s e n c eo fs i g n i fi c a n ts i d ee ff e c t s (Fig. 10-(b)).

Basis enrichment

As noticed above, the heat source reconstruction procedure is very shapesensitive. Indeed, Fig. 6-(c) shows that the reconstruction improves as the shape of the heat source gets closer to the mode shapes (e.g. curves rather than sharp angles). The basis enrichement procedure presented in [START_REF] Doudard | A probabilistic twoscale model for high cycle fatigue life predictions[END_REF][START_REF] Doudard | Identification of heat source fields from infrared thermography: Determination of 'self-heating' in a dual-phase steel by using a dog bone sample[END_REF] proposes to insert tailor-made modes into the natural basis in order to describe specifically the boundary conditions (and thus solve the aforementioned side effects problem). Yet, trivial mathematical considerations in finite dimension show that any orthogonal basis B which is added modes Qi looses its orthogonality property but still defines a unique projection of θ k (x, y)w i t h i n B = B S { Qi }.I n d e e d ,i ti ss h o w ni n [START_REF] Le Niliot | Infrared thermography applied to the resolution of inverse heat conduction problems: recovery of heat line sources and boundary conditions[END_REF]t h a tE q .( 7)i st h e solution that minimize the quadratic gap between the measured field and the projected approximation.

The hereby proposed procedure relies on the use of tailor-made modes (called thermomechanical modes) to describe the heat sources shapes themselves. A similar procedure has already been used in shape metrology with significant improvements [START_REF] Adragna | Analysis of shape deviations of measured geometries with a modal basis[END_REF]. The main difficulty is that users may only have an approximated knowledge of the heat sources shapes (for example from strain localization, geometrical specificities, etc) but the measured and projected data are temperature fields. Moreover, though the shapes of the heat sources can be roughly known a priori,t h e i rm a g n i t u d e sr e q u i r e st h e assumption of constitutive equations. Hence, the construction of thermomechanical mode requires (i) the use of a heat diffusion solver and (ii) a normalized heat mode. One can then use a Gram-Schmidt algorithm to orthogonalize the so-obtained base if needed but this is not mandatory since a projection within the dual basis (as proposed in Eq. ( 7)) is still available. In the presented numerical example, the addition of tailor-made modes that describes the heat sources in region no.1 and no.2 is performed from the actual knowledge of their shapes. Two enhanced bases are assembled using thermomechanical modes build from i) the exact heat sources shapes ( Q1 and Q2 )a n df r o mii) an approximation of it ( Qa 1 and Qa 2 ). Hence, the improvements achieved from using these enriched bases can be verified and quantified. In this section, only the results from enhanced modal projections are presented.

As shown in Fig. 12-(ab), four modes are built from the shapes of the heat sources no.1 and no.2 (two from the exact shape and two from the approximated shapes). Unitary heat sources are input and their temperature diffusion after 10 s are monitored (Fig. 12-(cd)). These modes are slipped after the first rigid mode (Q

1 =18(x, y) 2 R 2 )s ot h a t : BD = ⇣ Q 1 , Q1 , Q2 ,Q 2 ,...,Q N ⌘ Ba D = ⇣ Q 1 , Qa 1 , Qa 2 ,Q 2 ,...,Q N ⌘ . (26)
The so-obtained bases are reset orthogonaly using a Gram-Schmidt algorithm. Once again the temperature modes are normalized using an infinite norm. The approximated set of thermomechanical modes is built from eroded shapes of the heat sources. A Gaussian filtering-window is used over the exact heat source shapes (Fig. 12-(b)). This accounts for the fact that a priori knowledge of the heat source shapes may only be approximated. Fig. 13 shows the sum of reconstruction error versus signal-to-noise ratio. It can be seen that basis enrichment does improve the quality of the heat source calculation but only if the SNR is over 25 for BD and 32 for Ba D . Indeed, the natural basis remains the most able to describe the shape when measurement noise is massive. Another interesting feature depicted in Fig. 13 is the prevailing of exact modes Qi over the approximated ones Qa

i . An exact knowledge of the heat sources shapes leads to a closer assessment of their magnitudes. However, the modes built from an approximated knowledge of the shapes allows a better reconstruction than using the non-modified natural basis, in the case of small noise magnitudes. It can therefore be pointed out that, if the SNR is reasonnably high enough, any a priori knowledge of the heat sources shapes (approximated shapes and unitary magnitude) improves the heat source reconstruction. This can be seen and interpreted in Fig. 14.T h e m o d a l s p e c t r u m a r e depicted from projections at t =1 0 0s using B D , BD and Ba D .I tc a nb e seen that using an enriched basis leads to gather the high magnitude modes at very low frequencies (i.e. the left hand side of the spectrum). Accordingly the double space-derivation required by the Laplace term evaluation is not impaired by high frequencies modes and the numerical resolution provides a better match. As a matter of fact, the diffusion nature of the heat equation (1)a llows to build suitable temperature modes from rough estimate of the heat sources shapes (see Fig. 12-(cd)). In addition, one may think of plenty of formulations to enhance such an inverse thermo-mechanical problem espe-cially a combined approach using modes built from a priori knowledge of the heat sources shapes and/or boundary conditions, and/or noise distribution, and/or microstructural consideration, etc.I n d e e d ,t h eu s eo fc o u p l e d thermal-kinematic measurement offers a wide range of information on material deformation processes that could be used to enhance the projection procedure.

Conclusion

The presented study shows that regardless of how far from physics the modal projection is, it allows significant improvements in heat source estimation. It appears to be less sensitive to noise and side effects and thus allows a more accurate characterisation of calorific effects from temperature full-field measurements. In addition, the paper proposes to enhanced the projection bases by the means of thermomechanical modes possibly coming from strain measurements. This latter aspect has to be detailled more deeply in further publications. It is also noticable that building the modes from structual dynamics allows to account various geometry such as disk, cylinder, hollowed plates, notched specimen, etc since a very simple FE calculation is required.

Boundary conditions constitutes a classical problem of thermal field projection methods. It is noticable that the presented approach allows to impose awiderangeofthem. Indeedthemechanicalformulationoftheproblemleads to generate the modal basis from Finite Element calculations and hence allows to impose either displacement, forces or velocity which all have thermal equivalences.

The spectral basis, is built from thermal considerations and thus is closer to physics. However, it appears to be less efficient in the presented numerical tests. The formulation of this basis tested in this paper is the classical one and not a enhanced one such as proposed by [START_REF] Doudard | A probabilistic twoscale model for high cycle fatigue life predictions[END_REF][START_REF] Doudard | Identification of heat source fields from infrared thermography: Determination of 'self-heating' in a dual-phase steel by using a dog bone sample[END_REF]. Indeed, it should be accounted that several flaws of such a decomposition have already be pointed out in many former studies and improvement procedures has been proposed and implemented [START_REF] Neveu | Simulation de la conduction non linéaire en régime variable: décomposition sur les modes de branche[END_REF].

Though, such a complex filtering approach is probably unnecessary when measurements are clear from noise and exhibit sufficient space resolution, it constitutes a robust alternative in ill-posed inverse thermal problems. It must also be pointed that the presented numerical implementation leads to quick calculation since the modal basis is inverted once and for all. Processing a dozen images require less than a minute with an up-to-date computed.

Applying the presented procedure to actual experimental measurements is the main perspective of the study. Indeed, the presented investigation validates the ability of modal filtering approach to provides a robust and reliable estimation of the thermal heat sources.

  n e da b o v e ,w i t hb o u n d a r yc o n d i t i o n sp r e s c r i b e d on the four edges of the rectangle and no heat source. The thermal problem is solved, using Neumann conditions over ∂U,a s :

Figure 1 :

 1 Figure 1: Shape comparison between modal (a) and spectral (b) projection bases. Only the 25 first modes are presented.

Figure 2 :

 2 Figure 2: Solving flowchart highlighting the different steps of the proposed validation example.

Figure 3 :

 3 Figure 3: Meshes used for modes calculation and interpolation.

8 RMS 6 8 RMS 10 =7, 21 .10 8 RMS 2 8 RMS 6 8 RMS 10 =5

 86810218286810 =7,99.10 =2,33.10 =6,98.10

Figure 4 :

 4 Figure 4: Imposed (a) and computed (b-c) heat sources at various step time and related RMS errors (in W • m -3 ).

Figure 5 :

 5 Figure 5: Time evalution of the imposed and reconstructed heat sources for both projection bases.

Figure

  Figure6: Heat sources reconstruction using modal and spectral projection (red lines) for increasing noises compared with the imposed sources (black lines). Graphs depicte the evolution of the reconstructed heat sources along two lines that cross the main heat sources (region no.1 and no.2 presented in Fig.2).

Figure 7 :

 7 Figure 7: Comparison of the reconstruction error for increasing values of the signal-tonoise ratio.

Figure 8 :

 8 Figure 8: (a) Reconstruction error evolution versus the number of modes used for projection N . Four levels of noise are plotted. (b) Evolution of the optimum number of modes used in the projection versus signal-to-noise ratio.

Fig. 9 -Figure 9 :

 99 Figure 9: (a) Sketch of the FE model showing convective exchanges over the main surface. (b) Evolution of the RMS error versus h z for noiseless and noisy measurements.

Figure 10 :

 10 Figure 10: (a) Sketch of the FE model showing both heat sources and the heat tranfers over two sides of the domain U . (b) Shape of the x-derivative of the projected temperature fields in B ∆ and B D emphasizing side effects at t = 200 s.

Figure 11 :

 11 Figure 11: Evolution of the identified heat transfer coefficents (h x and h y ) versus signalto-noise ratio.

Figure 12 :

 12 Figure 12: (a) Exact source shapes. (b) Approximated source shapes. (c) Modes Q1 and Q2 corresponding to the temperature fields after 10 s diffusion. (d) Modes Qa 1 and Qa 2 corresponding to the temperature fields after 10 s diffusion of the approximated heat sources.

22

 22

Figure 13 :

 13 Figure 13: Sum of reconstruction error versus signal-to-noise ratio for B D , BD and Ba D .

Figure 14 :

 14 Figure 14: Projection spectrums at t = 100 s using natural basis and the two enhanced basis ( BD and Ba D ).


Selma Kchir

Tewfik Ziadi

Mikal Ziane

Serge Stinckwich
email: serge.stinckwich@ird.fr

A Top-Down Approach to Managing Variability in Robotics Algorithms

One of the defining features of the field of robotics is its breadth and heterogeneity. Unfortunately, despite the availability of several robotics middleware services, robotics software still fails to smoothly handle at least two kinds of variability: algorithmic variability and lower-level variability. The consequence is that implementations of algorithms are hard to understand and impacted by changes to lower-level details such as the choice or configuration of sensors or actuators. Moreover, when several algorithms or algorithmic variants are available it is difficult to compare and combine them.

In order to alleviate these problems we propose a top-down approach to express and implement robotics algorithms and families of algorithms so that they are both less dependent on lower-level details and easier to understand and combine. This approach goes top-down from the algorithms and shields them from lower-level details by introducing very high level abstractions atop the intermediate abstractions of robotics middleware. This approach is illustrated on 7 variants of the Bug family that were implemented using both laser and infrared sensors.

I. INTRODUCTION

The development of robotics software must deal with a large amount of variability from at least two sources. First, robots are very different from each other: "They have different locomotion mechanisms, different onboard computational hardware, different sensor systems, and different sizes and shapes." [START_REF] Smart | Is a Common Middleware for Robotics Possible?[END_REF]. Second, robots are used for very different tasks with varying constraints which leads to a large variety of algorithms.

Robotics middleware, among other improvements, has been a significant attempt at addressing the first kind of variability by decoupling robotics application from lower levels details. "It is designed to manage the heterogeneity of the hardware, [...] simplify software design [...]. A developer needs only to build the logic or algorithm as a component" [START_REF] Elkady | Robotics Middleware: A Comprehensive Literature Survey and Attribute-Based Bibliography[END_REF].

The task is huge however, as well explained by W.D. Smart who considers the hypothetical case of middleware providing an obstacle-avoidance routine for a mobile robot [START_REF] Smart | Is a Common Middleware for Robotics Possible?[END_REF]. According to him "writing a generic obstacle avoider that will work for all locomotion mechanisms, using input from all possible sensors is a daunting task". It is not surprising then that, a few years later, middleware services are still far from solving this decoupling problem.

Consequently, robotics algorithms are still

• difficult to understand,

• difficult to adapt or combine,

• impacted by changes in lower-level details.

Even with the very simplified assumptions such as those of the Bug family of navigation algorithms [START_REF] Lumelsky | Effect of Uncertainty on Continuous Path Planning for an Autonomous Vehicle[END_REF] it is far from obvious in which case such or such variant is best suited to go from one point to another while avoiding obstacles [START_REF] Ng | Performance Comparison of Bug Navigation Algorithms[END_REF]. In this paper we propose a top-down approach to complement the bottom-up middleware approach. The input of this approach is a robotic task and either a family of algorithms or at least enough knowledge to produce algorithms to solve it. Its output is twofold:

• a set of algorithmic, sensory or action abstractions, • a configurable generic algorithm. The produced algorithm can be configured by providing the values of a series of parameters to be adapted to different hypotheses, say on the environment. Furthermore it is generic in two other ways. First, it is decoupled from low-level details on sensors and actuators and second, the algorithmic variability which cannot be resolved statically by specifying configuration parameters is managed by dynamically linking the actions abstractions to executable routines. This rest of this paper is organized as follows. Section II describes our approach and illustrates it on Bug algorithms. Section III deals with the problem of organizing reusable implementations of the abstractions. Section IV discusses the approach. Section V presents preliminary validation results while section VI compares our approach to related work.

II. RATIONAL

In this section, we propose our approach to provide a generic algorithm for a family of algorithms for mobile robots. All the algorithms of the family are assumed to achieve the same task. It is further assumed that there is a significant plan aspect in the algorithms: they are not merely reactive. The generic algorithm can then be seen as a planner specialized to the considered task.

Our approach is to be performed by a human robotics expert with a strong background in programming. The input to the approach is a description of the task to be accomplished and a series of algorithmic descriptions in whatever format the expert understands (code, pseudo-code ...). The approach consists in defining a generic algorithm while extracting two kinds of abstractions to shield it from respectively low-level and algorithmic sources of variability. These two steps are performed concurrently. In a third subsequent step, the generic algorithm is implemented using the Template Method design pattern which delegates (similarly to the Bridge design pattern) sensing and actuation to so-called virtual sensors and actuators whose implementations can then vary independently from the algorithmic variants.

Figure 1 depicts the main steps of our approach (the first two steps are preformed concurrently):

1) Abstractions Identification

• Definition: Abstractions are a set of methods which allows to encapsulate specific data or to abstract specific methods by extracting a general signature from specific ones.

2) Generic Algorithm Definition

• Definition: A generic algorithm is a sequence of instructions written in terms of hardware and algorithmic abstractions.

• Decription: At this level, the robotic expert combine the abstractions identified previsouly in order to write the algorithm or to design a state machine corresponding to a specific robotic task. In this way, there is no specific detail related to a particular algorithms variant and the defined algorithm is completely independent from low level details and from specific algorithms variants.

3) Organising Implementation

• Definition: Organising implementation means implementing variants of a family of algorithms starting from the generic algorithm.

• Description: In this step, we implement the invariant parts of the algorithm in terms of abstractions. To do so, the Template Method (TM) design pattern [START_REF] Gamma | Design Patterns: Elements of reusable object-oriented software[END_REF] allows the definition of an algorithm skeleton in an abstract class as a template method. This method is the generic algorithm identified previsouly. It defines the basic step of the algorithm as "placeholders" methods (or hook methods), that are different in each subclass. Invariants parts of the algorithm are represented once in the abstract class as concrete methods. Abstract methods will be used to represent needed operations that are different in each subclass. The main idea is to represent the variants of the algorithm as subclasses that implement these abstract methods. Thus, algorithmic abstractions will be implemented in these subclasses. Concerning hardware abstractions, they are not implemented in subclasses like algorithmic abstractions but delegated to adaptors which implement them. We define one adapter per physical sensor. Adaptors implement these operations to extract physical data and convert them to the needed abstractions. Then, the specification of the physical sensor used in the algorithm is done at the deployment level. In the next section, we illustrate our approach on the Bug family of navigation algorithms.

III. CASE STUDY: BUG ALGORITHMS

It is symptomatic that even with the very simplified assumptions underlying these algorithms it is still difficult:

• to understand the differences among the algorithmic variants,

• to choose one variant over another,

• to share development efforts among variants.

A. Bug algorithms overview

The Bug algorithms attempt to solve the navigation problem in an unknown two-dimensional environment with fixed obstacles and a known goal position. Over 19 versions of Bug algorithms have been defined in the literature. Among them, 7 variants are considered in this paper.

Bug algorithms share the following assumptions [START_REF] Lumelsky | Effect of Uncertainty on Continuous Path Planning for an Autonomous Vehicle[END_REF] : 1) The environment is unknown and a finite number of fixed obstacles are placed arbitrarily. 2) The robot is considered as a point (i.e. without body). It has perfect sensors (for obstacle detection) and a perfect localization ability (e.g. to compute its distance to its goal).

In Bug family, the robot has a local knowledge and a global goal. In other words, the inputs of Bug algorithms are the robot's start position and the target position. At the end of the algorithm's execution, the robot must indicate if its goal is reached or if the goal is unreachable. Most of the Bug algorithms can be programmed on any mobile robot using tactile or distance sensing and a localization method while some require distance sensing.

B. Abstractions identification

In the litterature, as far as we know, Bug algorithms are only published with an informal description. Bug algorithms are very similar fundamentally but differ in some points of interest. Their principle is the following: [START_REF] Smart | Is a Common Middleware for Robotics Possible?[END_REF] The robot motion to its goal until an obstacle is detected on its way. (2) From the point where the obstacle were encountered (called hit point), the robot looks for a point (called leave point) around the encountered obstacle to be able to move to its goal again. (3) When a leave point is identified, the robot moves to it and leaves the obstacle. Steps (1), ( 2) and ( 3) are repeated until the goal is reached or until the robot indicates that the goal is unreachable. This description of Bug hides hardware and algorithmic variability. Hardware variability impacts on obstacle detection and localisation tasks because of their dependency of the robot actual sensors. Algorithmic variability is related to the leave point identification task. Following a goal directed approach and based on a comprehensive study of Bug algorithms, we have identified the following abstractions (described here as methods):

• Motion to goal: the robot needs to face its goal and to be able to go ahead it: FACEGOAL(POINT GOALPOSI-TION), GOAHEAD().

• Obstacle on the robot's way: BOOL OBSTACLEIN-FRONTOFTHEROBOT() -The local minimum condition. Using its distance sensors, the robot can detect discontinuity points [START_REF] Kamon | TangentBug: A Range-Sensor-Based Navigation Algorithm[END_REF] on obstacles on its way with respect to the target. -The disabling segment condition. A disabling segment occurs when the robot cannot move to its goal from all points in a segment while perfoming boundary following. -The step method. The robot uses its distance sensors to detect a point which is STEP [START_REF] Kamon | Sensory Based Motion Planning with Global Proofs[END_REF] closer to the target than any point already visited. Thus, we have identified the method FINDLEAVE-POINT(POINT ROBOTPOSITION, POINT HITPOINT, POINT GOALPOSITION) which will be applied to each point around the obstacle and which is specific to each algorithm variant to hide all low level details. The code of IDENTIFYLEAVEPOINT is given by the algorithm 1.

•
• Leave point identified: BOOL ISLEAVEPOINTFOUND()

C. Generic Bug algorithm definition

The generic algorithm is a combination of the previously defined abstractions as a sequence of instructions. Our generic algorithm is given in 2.

Algorithm 2: Bug generic algorithm

Sensors

: A perfect localization method.

D. Implementation: Template Method design pattern

As we said previously, the Template Method deals with variability problem by proposing to define the basic step of the algorithm as a template method written in terms of abstractions which will be implemented in subclasses. Define a subclass for each algorithm without handling sensors and actuators variability cause a combinatorial explosion. Consequently, we dealed with this problem by delegating hardware variability to what we call virtual sensors or adaptors. Virtual sensors convert specific data from the physical sensors to the needed data in the algorithm. They implement a set of interface abstract operations. For instance, the method OB-STACLEINFRONTOFTHEROBOT() is different in each sensor adaptor. We have defined an adaptor per sensor. In case of multiple physical sensors use, their adaptors are combined together to provide the needed data for the algorithm. The architecture of our implementation is presented in Figure 2.

To perform wall following, right hand and left hand algorithms are written in terms of sensors abstractions (i.e. OBSTACLEONTHELEFT() and OBSTACLEONTHERIGHT())

E. Discussion

The hierarchy of classes defined by our approach can lead to a combinatorial explosion if we add additional variants of Bug family. It is then much cheaper to build operations of the algorithm as components and to assemble desired family members from them. This is our main motivation for using an alternate solution based on Software Product Lines (SPL). SPL are easy to use, compact and generate only a configuration which interests the user. This constitutes an immediate topic for a future work.

IV. RESULTS AND VALIDATION

Several performance comparison studies [START_REF] Ng | Performance Comparison of Bug Navigation Algorithms[END_REF] were realized on the Bug family. In this section, we do not intend to perform a comparison between Bug family but to prove that each algorithm of the Bug family fits with our generic algorithm. We demonstrate the capabilities of our generic algorithm in the OROCOS-RTT framework through 7 variants of Bug: Bug1 [START_REF] Lumelsky | Effect of Uncertainty on Continuous Path Planning for an Autonomous Vehicle[END_REF], Bug2 [START_REF] Lumelsky | Dynamic Path Planning for a Mobile Automaton with Limited Information on the Environment[END_REF], Alg1 [START_REF] Noborio | A Comparative Study of Sensor-based Path-planning Algorithms in an Unknown Maze[END_REF], Alg2 [START_REF] Sankaranarayanar | Path Planning for Moving a Point Object amidst Unknown Obstacles in a Plane: a New Algorithm and a General Theory for Algorithm Development[END_REF], Dist-Bug [START_REF] Kamon | Sensory Based Motion Planning with Global Proofs[END_REF],TangentBug [START_REF] Kamon | TangentBug: A Range-Sensor-Based Navigation Algorithm[END_REF] and Rev1 [START_REF] Horiuchi | Evaluation of Path Length Made in Sensor-Based Path-Planning with the Alternative Following[END_REF].

Simulation was performed using Stage-ROS with 2 configurations. The first tested configuration was done with a laser scanner with 180 degrees scanning angle and a detection range which varies from 0.02 meters to approximately 4 meters and a GPS for localization. The robot does not have any knowledge about its environment except its start position and the goal position. The second configuration was done using 3 infrared range sensors placed on the front, on the left and on the right compared to the central axis of the robot with a field of view equals to 26 degrees and a range which detects until 2 meters. To demonstrate that any algorithm of the ones studied here fits with our generic algorithm, we have tested both configurations in 3 different environments for all algorithms. We defined the first environment (see figure 3) to validate the target reachability condition. In all implemented algorithms, the robot returns failure because it can not achieve its goal. The second environment shown in figure 4 is a simple simulation environment with one obstacle. Some algorithms behave similarly in this environment despite their different obstacle avoidance strategies. For instance, Alg1 and Bug2 rely on the M-line detection condition but behave differently when they encounter an already traversed point. In other words, Alg1 was defined above Bug2 to overcome situations where the robot can find itself in an overlasting loop around obstacle. For this reason, we defined the third environment, presented in figure 5, to validate the encountered points condition (i.e. when the robot encounters an already traversed point), particularly used in the algorithms Alg1 and Alg2. The trajectory of the robot could depend on the computation time of sensors information. Since most of the algorithms treatments are reactive, we had to define the period of execution time which allow us to get as much realtime information as possible and to optimize the execution time and the robot's path. Consequently, we set the execution period of the algorithm to 0.5 seconds.

Alll the code of Bugs algorithm variants1 and simulations results2 are available on GitHub.

V. RELATED WORK Several software engineering technologies and methods aim at improving software design and reusability. In robotics, reusability is obtained after handling variability which could be related to the robot hardware or to the robot's capabilities (e.g. different strategies for obstacle avoidance, etc). Robotics software like Robot Operating System (ROS) [START_REF] Quigley | ROS: an Open-source Robot Operating System[END_REF], RISCWare framework [START_REF] Elkady | A Structured Approach for Modular Design in Robotics and Automation Environments[END_REF], Player [START_REF] Collett | Player 2.0: Toward a Practical Robot Programming Framework[END_REF] propose a hardware abstraction layer to encapsulate the sensors that gather information about the environment and provide a set of predefined intefaces. These middlewares rely on a bottom up approach which consists on classifying the most used physical sensors (or actuators), analysing the potential data they are able to provide and then define interfaces. Unlike these middlewares, we use a top-down approach which consists in analysing the needed data in a particular application, defining abstractions and then writing the application in terms of these abstractions (possibly provided by middlewares).

In our approach, we have applied the Template Method design pattern to implement variants of a generic algorithm. They are several proposals that try to apply design patterns in robotics to handle variabilities. CLARATy authors [START_REF] Nesnas | CLARAty: Challenges and Steps Toward Reusable Robotic Software[END_REF] say they used many well-know techniques developed by the software community including design patterns but without being much more explicit about them. MARIE [START_REF] Côté | Robotic Software Integration Using MARIE[END_REF], a middleware framework for robotics, applied the Mediator Design Pattern [START_REF] Gamma | Design Patterns: Elements of reusable object-oriented software[END_REF] to create a mediator interoperability layer for distributed robotics applications.

VI. CONCLUSION AND PERSPECTIVES

In this paper we have proposed an approach to organize families of algorithms so that the algorithmic decisions to make are clearly expressed and decoupled from implementation details. Our approach relies on the Template Method design pattern to define a generic algorithm for Bug algorithms family.

Our approach consists of three steps. The first step takes as input a set of algorithmic variants; as described in the litterature, and manually extract hardware abstractions and what we have called algorithmic abstractions. The generic algorithm is then defined as a sequence of instructions in terms of these abstractions. The implementation of these variants relies then on the Template Method design pattern.

This approach has been illustrated on the Bug family of robot navigation algorithms. Seven implemented variants of Bug algorithms have been implemented using the OROCOS-RTT robotic framework. Simulation was performed in different unknown environments with a random positioning of obstacles.

Bug variants are about 20 versions of algorithms. Implementing all BugAlgorithm subclasses is very complicated to understand and too expensive to build all Bug family members. In addition, it is error prone to not define constraints on subclasses and the data they must specify and use. It is then much cheaper to build operations of the algorithm as components and to assemble desired family members from them. This is our main motivation for using Software Product Lines (SPL) [START_REF] Clements | Software Product Lines: Practices and Patterns[END_REF] as an immediate topic of future work.

Fig. 1 .

 1 Fig. 1. Generic algorithm: from identification to implementation

  Current position: POINT GETPOSITION() • Get around the encountered obstacle consists in performing clockwise or anticlockwise circumnavigation while maintaining a safety distance from the obstacle: WALLFOLLOWING(DIRECTION) -DOUBLE GETSAFEDISTANCE() -BOOL OBSTACLEONTHERIGHT(), BOOL OBSTA-CLEONTHELEFT(): Depending if we follow the wall on the right or the left. -DOUBLE GETRIGHTDISTANCE(), DOUBLE GETLEFTDISTANCE(): To keep a safety distance from a wall, we need to know the distance to the wall from the right and the distance from the left. • Look for a leave point: IDENTIFYLEAVEPOINT(BOOL DIRECTION, POINT ROBOTPOSITION, POINT GOAL-POSITION) which consists in wall following while looking for a leave point. This could be a local decision (choose the first leave point which satisfies the algorithm condition) or a global decision (choose the best leave point after visiting all the points around the obstacle). The following conditions are examples of decisions that the algorithm defines to find a leave point: -The closest point around obstacle boundary condition. It consists on recording the closest point to the goal among all the points ever visited by the robot while performing boundary following[3][6]. -The m-line detection condition. It is a straight line between the starting point and the target which aims at providing a set of prededefined leave points. The robot must leave only on these points[6][7].

Fig. 2 .

 2 Fig. 2. Bug algorithms class diagram (extract)

Fig. 3 .Fig. 5 .

 35 Fig. 3. Environment with unreachable target Fig. 4. Basic Simulation Environment

https://github.com/SelmaKchir/BugAlgorithms

https://github.com/SelmaKchir/BugAlgorithms/wiki/Implementing-Bug-Algorithms-variants


est destinée au dépôt

Palaeoenvironments and palaeoclimate records in lacustrine deposits at the Eocene-Oligocene Transition. New insight from de CDB1 borehole (Rennes Basin, France)

Hugues Bauer, Jean-Jacques Châteauneuf, Chantal Bourdillon, Pierre

Saint-Marc, François Guillocheau

To cite this version:

Hugues Bauer, Jean-Jacques Châteauneuf, Chantal Bourdillon, Pierre Saint-Marc, François Guillocheau. Palaeoenvironments and palaeoclimate records in lacustrine deposits at the Eocene-Oligocene Transition. New insight from de CDB1 borehole (Rennes Basin, France

00995163 Palaeoenvironments and palaeoclimate records in lacustrine deposits at the Eocene- Oligocene Transition.New insight from de CDB1 borehole (Rennes Basin, France)

  The Chartres- corresponds to alternating sandy and clayey deposits. Depositional environments range from fluvio-lacustrine to fluvio-marine settings with occasional mangal development as attested by Avicennia.The formation is assigned to the Bartonian by benthic foraminifera and palynology. The lowermost samples yielded a "Biarritzian" age, which is equivalent to the Early-Middle Bartonian. The palynological record shows a gradual palaeoclimate change. The Bartonian is quite similar to the Lutetian of the Paris Basin, with a warm and humid "tropical" climate. The Early Priabonian, up to 278m, is still under humid and warm conditions, but development of herbaceous vegetation attests for a slight seasonality. During the Late Priabonian, a large development of Pinaceae coeval with a decrease in megathermic flora points to the EOT major climatic change. Indeed, Early Rupelian assemblages show even greater percentages in Pinaceae (above 50%) despite a steady lacustrine environment. The very last megathermic elements disappear at the base of the Natica Marls Fm. The clay mineralogy evolution looks more abrupt. The clay assemblage from the Bartonian to the Late Priabonian is fully kaolinitic. The Early Rupelian assemblage is still dominated by the kaolinite but smectite appears in a significant amount (up to 60%). Whether the sudden

	Bauer H. 1 , Châteauneuf J.J.C. 2 , Bourdillon C. 3 , Saint-Marc P. 4 &Guillocheau F. 5
	1 BRGM, ISTO, UMR 7327, BP 36009, 45060 Orléans, France	h.bauer@brgm.fr
	2 BRGM (retired) 8 quai du Châtelet, 45000 Orléans, France	chateauboy@wanadoo.fr
	3 ERADATA, 170 avenue FélixGeneslay 72100 Le Mans, France bourdillon@eradata.fr
	4 SEMM LOGGING, Les Maufras -18360 Vesdun, France	saint-marc@semmlogging.fr
	5 Université de Rennes 1, UMR 6118 Géosciences Rennes, Campus de Beaulieu, 35042
	Rennes Cedex, France	francois.guillocheau@univ-rennes1.fr
	On behalf the CINERGY project, a 675m-long borehole was cored through the sedimentary
	column (405m) and the basement (270m).
	The studied interval encompasses 3 lithostratigraphic units, from 66m to 405m depth. The
	chronostratigraphic framework of the sedimentary series relies on benthic foraminifera and
	palynology. Palaeoenvironmental reconstruction is based on sedimentology, pollen analysis
	and clay mineralogy.

). 19th International Sedimentological Congress : ISC 2014, Aug 2014, Genève, Switzerland. hal-The Natica Marls Fm (66-85m) is a lagoon-marine unit exhibiting metric sequences from restricted bay to salt marsh (schorre). Pollen grains, Mollusks and Dinocysts assemblages give a Rupelian age. The Lower SapropelsFm (85-375m) is the thickest unit and exclusively made of lacustrine and palustrine clay deposits. Clays are either thinly laminated (varve-like) or massive, blocky with pedoturbation and/or brecciated fabric, alternating in thick (20-40m) sedimentary sequences. Both facies show varying organic content, up to 40% TOC. Pollen assemblages show a bimodal repartition between the laminated and massive facies. The former are interpreted as an open lacustrine system and the latter, as a closed lacustrine system whose floating mats vegetation, characterized by papyrus and lotus is typical of permanent flooded areas. The first occurrence of the Early Rupelian marker B. hohli is observed at 195.13m. First results from magnetostratigraphy and cyclostratigraphy argue for an E/O boundary at ~202 m depth.The Eocene-Oligocene Transition is thus recorded in a detailed, continuous depositional environment.

Acknowledgements: This work was funded by Région Bretagne, CG35, Chartres-de-Bretagne, Rennes Métropole, IAV, SMG35, SMPBR, Préfecture de Bretagne, the BRGM, AELB and ADEME agencies.


Hugues Bauer

Robert Wyns

Paleoweathering profile developed on homogenous sedimentary basement: an integrated approach from the CDB1 deep borehole (Rennes Basin, France)

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Paleoweathering profile developed on homogenous sedimentary basement: an integrated approach from the CDB1 deep borehole (Rennes Basin, France).

Bauer H. 1 , Wyns R. 1 1 BRGM, ISTO, UMR 7327, BP 36009, 45060 Orléans, France. h.bauer@brgm.frr.wyns@brgm.fr Weathering profiles are mostly studied on their upper part (crust, saprolite) where leaching and concentration/precipitation of valuable element occur. For water resource and hydrocarbon purpose, the transition between saprolite and fresh basement is of utmost importance. Here is found the fissured layer, a highly fissured bedrock that it favorable for aquifer/reservoir properties. Such aweathering profile model is well known on igneous and metamorphic rocks, but curiously very poorly documented on sedimentary, non-metamorphic, basement. On behalf the CINERGY project, a 675m-long borehole (CDB1) was cored through the sedimentary infill (405 m), then the weathered basement (110 m) and finally the fresh bedrock (160 m). The basement is made of the "brioverianschists", which actually are epimetamorphicshales and fine grained sandstones, highly folded and cleaved.

Here is presented an integrated (well-logging, mineralogy, petrography) study of the weathering profile that has been preserved under the sedimentary infill of the Rennes Basin. The well logging tools include Gamma-Ray, long and short resistivity, neutron porosity, gamma-gamma density, Pef, Full wave Sonic, BHTV and caliper. The cores were regularly sampled for total rock and clay mineralogy and some petrographic control were realized on specific facies and fracture fillings. The core description gives a first visual sequence of the weathering profile. From up to bottom, we observe a 1 m-thick interval of massive ochrish clays, followed by 0.7 m-thick of structured ochrish clays, becoming greyish downward. The clays tend to become harder downward, and from around 422 m to ~470 m, the shales look highly fractured. From 470 m to the bottom of the core, the rocks look very hard. When looking at the physical parameters from the well logging, the interval boundaries do not look the same, and especially the saprolite/fissured layer boundary is very hard to point out, as the parameters show a gradual change all along the weathered profile. It is only towards 520 m that parameters (GR, sonic, resistivity) seem to reach a steady state, highlighted by the very monotonous lithology of the basement. The clay mineralogy is made of kaolinite, illite and chlorite. Chlorite is dominant below 510 m whereas kaolinite becomes more abundant from this depth upward; illite looks quite constant throughout the core. Petrographic control on the fracture fillings shows a polyphased precipitation of chloritic quartz, followed by carbonated (calcitic and dolomitic) and pyritized fillings. Oxidization of pyrite into goethite is observed even deep into the weathering profile (480 m). These first results show that the classic model of weathering profile is not suitable for sedimentary basement, which is in this case made of argillaceous sediments which did not develop a real fissured layer and so did not allow the aquifer/reservoir layer to take place.

Acknowledgements: The CDB1 borehole was funded by Région Bretagne, CG35, Chartresde-Bretagne, Rennes Métropole, IAV, SMG35, SMPBR and Préfecture de Bretagne. The BRGM, AELB and ADEME agencies funded the scientific survey and geological studies.