Rada Alhasan

Thierry Lubin

Jean Lévêque

Kévin Berger

Bruno Douine

Abderrezak Rezzoug

Smail Mezani

Gaëtan Didier

Melika Hinaje

Denis Netter

Study of a superconducting motor with high specific torque

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche

Introduction

Improving the performance of electrical machines is a big challenge for decades [START_REF] Tixador | 50-kW experimental superconducting permanent-magnet motor[END_REF][START_REF] Netter | Theoretical study of a new kind HTS motor[END_REF][START_REF] Kalsi | Development status of superconducting rotating machines[END_REF][START_REF] Masson | Experimental study of a new kind of superconducting inductor[END_REF][START_REF] Jo | Characteristic study and three dimensional magnetic field analysis of the superconducting synchronous machine[END_REF]. One of the most important parameters to improve machine performance is the generation of a high magnetic field by an inductor. Superconductors are very attractive materials because of their absence of resistive losses when they are fed by DC current, and due to the very high current densities. For these reasons, they have the advantage of allowing to build lighter electrical machines with better performance than conventional motors. We propose a new superconducting inductor structure. The topology uses an inductor formed by two superconducting field coils and a bulk superconductor which is a magnetic screen leading the flux density in the air gap . A three-dimensional analysis of the magnetic field, is made, using a 3D finite element software, to study the performance of the proposed inductor.

Principle of the superconducting inductor

Figure .2 shows the structure of the proposed inductor. It is composed of two coaxial superconducting coils. The coils are supplied with electric currents having the same direction to generate a magnetic field (B1, B2). A bulk superconductor is placed between the two coils. The superconductor is inclined along the length of the inductor. This superconductor is used as a magnetic barrier, and is located between the two solenoids, thereby creating a variable magnetic field. In order to improve the radial distribution of induction, we also used a ferromagnetic material between the two coils. The superconductor is inserted between the ferromagnetic parts as shown in Figure 1, Iron guide the magnetic field toward the centre of the inductor and therefore reduce the leakage field near the coil.

The geometric parameters of the proposed structure are given in Table 1 These parameters were chosen after a sizing of a superconducting motor to obtain a high magnetic field and correspond to a realistic size machine. For the simulations, a current density of 650A / mm2 will be chosen by using the behaviour law for low critical superconducting temperature wire.

Study of the whole superconducting motor

The principle of the whole superconducting motor is described on figure 2. We have a classical cooper armature or a superconducting armature with a back iron. In the first case, we have one cryostat for the superconducting inductor, in the second case, we need to design a cryostat for the inductor and the armature. In this paper we present results with a classical copper armature and a superconducting inductor with an airgap of one centimetres. This airgap is very high due to the cryostat which are located inside. So, after simulation using FEM software, we obtain the distribution of magnetic flux presented on figure 3.

The magnetic field varies between -4T and 4T. If we compare to a classical motor it is four time greater. Therefore, we can increase by a factor 4 the specific power and torque of the motor. This magnetic field is calculated on the outer radius of the inductor

Experimental results

We decided to build a machine based on this principle. The superconducting coils will be made in NbTi. The NbTi superconducting wire is a low critical temperature material cooled by liquid helium at 4.2 K. The superconducting screen is an YBCO plate. Initially, before realizing the entire machine, we tested the inductor alone. To this end, all of the inductor, coil and YBCO bulk, are immersed in nitrogen tank. The YBCO bulk is superconducting at the temperature of liquid nitrogen. The coils made of NbTi, will not superconducting at this temperature. They will therefore be fed with a very low current to avoid destructive heating of the wire but still enough to have a measuring magnetic field. The purpose of this experiment is to check the form of magnetic induction created by our induction to verify our computational tools. The inductor is realized with geometric parameters given in the table 1, but with a current density of 10.6 A/mm2. The experimental bench consists of our inductor, and a Hall probe that we can move along the height of the inductor and along its bore diameter. We also have the supply of the superconducting coils, of the Hallprobe and the acquisition system. This is shown in Figure 4. The results are shown in the figures 5 to 7 and show the correlation value of shape of the magnetic field between the theoretical and experimental results. The measurements are made at a distance of 1.5 cm from the inductor. In Figure 6, we have shown the variation of the induction along the bore diameter. In Figure 7, we have shown the variation of the induction along the useful length of the inductor. The observed induction values are low because we are far from the current rating of the inductor, 4 A instead of 200 A, when the superconducting coils will be in liquid helium.

Fig. 1 :

 1 Fig. 1: Principle of the superconducting inductor

Fig. 2 :

 2 Fig. 2: Principle of the superconducting motor

Fig. 3 :

 3 Fig. 3: Magnetic flux in the airgap of the superconducting motor

Fig. 4 :

 4 Fig. 4: Experimental bench

Fig. 5 :Fig. 6 :Fig. 7 :

 567 Fig. 5: Shape of the magnetic field

Table 1 :

 1 Geometrical parameters of the inductor

	Outer radius	Rex (mm)	68
	Inner radius	Ri (mm)	53
	Distance between solenoids C (mm)	45
	Solenoid length	L (mm)	45
	Thickness of the screen	E (mm)	10
	Density de courant	J (A/mm 2 )	650


Bruno Després
email: bruno.despres@cea.fr

Raphaël Loubère
email: loubere@mip.ups-tlse.fr

Convergence and Sensitivity Analysis of Repair Algorithms in 1D

Keywords: Finite Volume Schemes, Repair algorithms, TVD schemes, Gas Dynamics

We prove the convergence of some repair algorithms for linear advection in dimension one. The convergence depends on the size of the box where the distribution of the mass excess is performed. Various numerical examples illustrate the theoretical results. Applications to gas dynamics in dimension one is also discussed.

Introduction

The aim of this work is the convergence analysis of some very simple Repair Algorithms. To our knowledge such a proof has never been given, even for simple problems. The Repair algorithms we consider are very much in the spirit of what was proposed by M. Shashkov and B. Wendroff in [START_REF] Shashkov | Wendroff -The Repair Paradigm and Application to Conservation Laws[END_REF]. The underlying idea of Repair methods is related to the fact that conservative remapping methods (present in Arbitrary-Lagrangian-Eulerian (ALE) framework for example) may not necessarily preserve a maximum principle: no new extrema should be generated during the remapping process. If so a Repair method can remove these new maxima (resp. minima) by distributing (resp. taking) the amount of the associated conservative variable to (resp. from) the neighborhood. In this paper we limit the theoretical study to transport equation in 1D and to some simple remapping strategy in 1D. However we hope that some of the conclusions of the study can give insights even for multidimensional repair algorithms for systems of PDE's, such as compressible gas dynamics. The outline of this paper is as follows. In the second section, we describe the repair algorithm for the linear transport. In section three we adapt the repair algorithm for the remapping process. The numerical results on 1D linear transport are gathered in the fourth section, whereas the fifth section presents some numerical results for 1D gas dynamics. Finally the conclusion and the plan for future developments is presented in the last section.

Repair algorithm for linear transport

Let us consider the equation of transport in 1D, with constant velocity ∂ t u + a∂ x u = 0, a > 0.

(

) 1
The initial condition is u(0, x) = u 0 (x) where u 0 ∈ L 1 (R) ∩ BV (R). Let us consider a mesh, that is uniform even if our result extends to non uniform mesh. The mesh size is denoted by ∆x. The time step is denoted by ∆t. Thus u n j stands for the numerical solution in cell j at time step n.

The repair strategy amounts to: 1) compute a prediction of the numerical solution at time step n + 1 using a "reasonable" and "local" scheme, as instance this scheme can be a high order non monotone scheme or a highly anti-dissipative scheme, 2) check if the new value satisfies a local maximum principle, 3) if the new value does not satisfy the local maximum principle, then repair it. Repairing means changing the value of the unknown for 3) to be fulfilled. A difficulty is that one wants the total mass to be preserved. So we need to describe in details how to redistribute the mass of the repaired quantity, such as the total mass is preserved.

This family of algorithms can be local if one redistributes the mass in a local box around the cell that needs to be repaired, or global if one redistributes the mass in the entire domain. It has been a debate since the early time of the Repair Algorithms to decide whether local repairing is better or not than global repairing. On one hand, global repairing is a more simple algorithm than local repairing. On the other hand our analysis implies that we can not prove the convergence of a global repair process. This is a theoretical indication that global repair process can be dangerous in some cases. Moreover a global repair process breaks the causality by instantaneously spreading mass all over the domain. The numerical experiments show it is indeed the case, in particular when the prediction step or the underlying physics is oscillating.

Let us describe in detail the repair algorithm that we analyze in this work.

The prediction scheme

First, one computes the new value of the unknown using the finite volume and conservative scheme

u n+1 j -u n j ∆t + a u n j + c n j+ 1 2 -u n j-1 -c n j-1 2 ∆x = 0, (2)
International Journal on Finite Volumes where u n j is the mean value of u(x) in cell j:

[x x-1 2 ; x x+ 1 2
]. The Courant or CFL number is less than one (a > 0)

ν = a ∆t ∆x ≤ 1. ( 3
)
The scheme is a priori different from the upwind scheme. Thus all the difference between the upwind scheme and the scheme used is embedded in the definition of the correction flux ac n j+ 1 2 for all j. This correction flux can be either a linear or a non linear function of (u n j ). The only assumption is

∃C 1 > 0, ∃k ∈ N, |c n j+ 1 2 | ≤ C 1 j-k≤q≤j+k |u q -u q-1 |. (4)
The hypothesis (4) essentially means that the flux is defined as the upwind flux plus a correction. Of course the correction is zero if the numerical profile is flat (that is if u qu q-1 ≡ 0 in a neighborhood of cell j) : ( 4) is compatible with such a principle. The hypothesis is true for the Lax-Wendroff scheme as instance, moreover all non linear TVD algorithms satisfy (4).

Correction step

The spirit of this repair algorithm is to compare u n+1 j with

M n j = max(u n j , u n j-1 ) and m n j = min(u n j , u n j-1 ), (5)
that is one checks if m n j ≤ u n+1 j ≤ M n j is true or not. Suppose u n+1 j > M n j , then one has to modify the value of u n+1 j and redistribute the mass "around". In the convergence analysis of the method, we discovered that it is better at the theoretical level not to redistribute the mass globally but locally at least in a box of size p ∈ N around the current cell. This is why we have introduced a new step in the repair algorithm first proposed in [START_REF] Shashkov | Wendroff -The Repair Paradigm and Application to Conservation Laws[END_REF] to be able to ensure that the redistribution of mass can be made in the box of size p. Since p is a parameter of the method, one recovers the global repair by setting p ≈ +∞.

So let us define boxes of size p. Each box is the collection of cells j such that rp ≤ j ≤ (r + 1)p -1 where r ∈ Z. The mathematical definition of these boxes B r is

B r = {j; rp ≤ j ≤ (r + 1)p -1}, r ∈ Z. ( 6
)
However it is also possible to use boxes of different sizes, provided the size is smaller than the predefined maximal box's size p. It is also possible to use moving boxes, that is the starting point of each box is different from one time step to the other. For the simplicity of the mathematical exposure we use only [START_REF] Després | Lagoutière -Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF].

To make the correction we first need to compute

b M r =   j∈Br u n j -ν(u n j -u n j-1 ) -M n j   -ν(c n (r+1)p-1 2 -c n rp-1 2 ), (7) b
m r =   j∈Br u n j -ν(u n j -u n j-1 ) -m n j   -ν(c n (r+1)p-1 2 -c n rp-1 2 ). ( 8
)
International Journal on Finite Volumes

Most presumably b M r ≤ 0 (resp. b m r ≥ 0), since this is the result of a comparison be-

tween j∈Br u n j -ν(u n j -u n j-1 ) -M n j ≤ 0 (resp. j∈Br u n j -ν(u n j -u n j-1 ) -m n j ≥ 0) and ν(c n (r+1)p-1 2 -c n rp-1 2
). Moreover if p is large enough and ν small enough (i.e. the time step is small) then b M r ≤ 0 (resp. b m r ≤ 0) is probably stating true. The correction is here to ensure that b M r ≤ 0 and b m r ≥ 0 are always satisfied. The idea being that if these inequalities are not satisfied then we multiply the value of the fluxes by a small number such that b M r ≤ 0 and b m r ≥ 0 are fulfilled. So we define

d n j-1 2 = µ n j-1 2 c n j-1 2 where µ n j-1 2 ∈ [0, 1]. (9)
The coefficient µ n j- is almost equal to the flux of the prediction scheme. We

need to check bM r -ν(µ n (r+1)p-1 2 c n (r+1)p-1 2 -µ n rp-1 2 c n rp-1 2 ) ≤ 0, ( 10
) bm r -ν(µ n (r+1)p-1 2 c n (r+1)p-1 2 -µ n rp-1 2 c n rp-1 2 ) ≥ 0, (11)
where by definition bM r = j∈Br u n jν(u n ju n j-1 ) -M n j is not positive and bm r = j∈Br u n jν(u n ju n j-1 )m n j is not negative. One feasible strategy can be derived as:

International Journal on Finite Volumes Analysis of (10) if c n (r+1)p-1 2 ≥ 0 and c n rp-1 2 ≥ 0 : Then (10) is true once ν(µ n rp-1 2 c n rp-1 2 ) ≤ -bM r . Thus we define ϕ 1,- r , ϕ 1,+ r such that µ n rp-1 2 ≤ ϕ 1,- r = -bM r c n rp-1 2 and ϕ 1,+ r = +∞. ( 12
) if c n (r+1)p-1 2 < 0 and c n rp-1 2 < 0 : Then (10) is true once ν(-µ n (r+1)p-1 2 c n (r+1)p-1 2
) ≤ -bM r . Thus we define ϕ 2,- r , ϕ 2,+ r such that ) ≤ bm r . Thus we impose that

ϕ 2,- r = +∞ and µ n (r+1)p-1 2 ≤ ϕ 2,+ r = -bM r -c n (r+1)p-1 2 . ( 13
) if c n (r+1)p-
ψ 1,- r = +∞ and µ n (r+1)p-1 2 ≤ ψ 1,+ r = bm r c n (r+1)p-1 2 . ( 15
) if c n (r+1)p-1 2 < 0 and c n rp-1 2 < 0 : Then (11) is true once ν(-µ n rp-1 2 c n rp- 1 2
) ≤ bm r . Thus we impose that

µ n rp-1 2 ≤ ψ 2,- r = bm r -c n rp-1 2 and ψ 2,+ r = +∞. ( 16
) if c n (r+1)p-1 2 ≥ 0 and c n rp-1 2
< 0 : Then it is not possible to simplify the inequality [START_REF] Shashkov | Wendroff -The Repair Paradigm and Application to Conservation Laws[END_REF]. Thus we impose ψ 3,- (18) Inside the boxes we do not correct, that is

r = ψ 3,+ r that µ n (r+1)p-1 2 , µ n rp-1 2 ≤ ψ 3,- r = -bM r -c n (r+1)p-1 2 + c n rp-1 2 (17) if c n (r+1)p-
d n j-1 2 = c n j-1 2 , ∀j = rp. ( 19
)
The next step consists in the computation of the new prediction ûn+1 j with the corrected flux:

ûn+1 j -u n j ∆t + a u n j + d n j+ 1 2 -u n j-1 -d n j-1 2 ∆x = 0. ( 20
)
International Journal on Finite Volumes Lemma 2.2 One has the inequalities after the correction step

  j∈Br u n j -ν(u n j -u n j-1 ) -M n j   -ν d n (r+1)p-1 2 -d n rp-1 2 ≤ 0 (21)   j∈Br u n j -ν(u n j -u n j-1 ) -m n j   -ν d n (r+1)p-1 2 -d n rp-1 2 ≥ 0. ( 22
)
The proof is performed by considering all cases (12-17) separately.

Inequalities (21-22) will be crutial in the analysis of the repairing procedure, essentialy inequality (31) in subsection 2.4.

Repairing

As already mentioned, to repair a value means successively to compare with the local maximum and minimum, to truncate if needed, then to redistribute the excess of mass on all surrounding cells. Using mathematical notations, one gets

   if ûn+1 j > M n j then u n+1 j = M n j and ∆m n j = ûn+1 j -M n j > 0, if ûn+1 j < m n j then u n+1 j = m n j and ∆m n j = ûn+1 j -m n j < 0, else then u n+1 j = ûn+1 j and ∆m n j = 0. (23)
The total mass in box B r of the new unknown u n+1 j may as well be different from the correct mass, so one defines the default of mass as:

∆M n r = rp≤j≤(r+1)p-1 ∆m n j . ( 24
)
This default of mass may be positive or negative. So we need to redistribute it on the box to get at least a conservative algorithm. Following [START_REF] Shashkov | Wendroff -The Repair Paradigm and Application to Conservation Laws[END_REF], we consider

   if ∆M n r > 0 then u n+1 j = u n+1 j + λ n r (M n j -u n+1 j ), if ∆M n r < 0 then u n+1 j = u n+1 j + λ n r (m n j -u n+1 j ), if ∆M n r = 0 then u n+1 j = u n+1 j , (25)
where the coefficient λ n r is set to

       if ∆M n r > 0 λ n r = ∆M n r P rp≤j≤(r+1)p-1 (M n j -u n+1 j ) , if ∆M n r < 0 λ n r = ∆M n r P rp≤j≤(r+1)p-1 (m n j -u n+1 j ) , if ∆M n r = 0 λ n r = 0. ( 26
)
The repair algorithm that we analyze in this paper consists of equations ( 2) to (26).

International Journal on Finite Volumes

Properties

Before proving our main convergence theorem, we state the stability lemma Lemma 2.3 Whatever the value of p ∈ N ⋆ is, for all time step n, the repair algorithm is such that the total mass is preserved (conservation)

j u n+1 j = j u n j , (27)
the maxima and minima are respected (maximum principle)

m n j ≤ u n+1 j ≤ M n j , ∀j, (28)
and the total variation is diminishing (TVD)

j |u n+1 j -u n+1 j-1 | ≤ j |u n j -u n j-1 |. ( 29
)
Remark: One notices that the scheme is TVD because of (29). Proof First one has

j∈Br ûn+1 j =   j∈Br (u n j -ν(u n j -u n j-1 ))   -ν d n (r+1)p-1 2 -d n rp-1 2 .
Since by construction u n+1 j = ûn+1 j -∆m n j , then we get the relation

∆M n r + j∈Br u n+1 j =   j∈Br (u n j -ν(u n j -u n j-1 ))   -ν d n (r+1)p-1 2 -d n rp-1 2 . ( 30
)
Assume for instance that ∆M n r > 0 (the other cases are easily deduced by mimicking this one). Then the next stage of the algorithm consists in the computation of λ n r . The key point is to prove that 0 ≤ λ n r ≤ 1 and the property follows. One has

∆M n r -j∈Br M n j -u n+1 j = j∈Br (u n j -ν(u n j -u n j-1 ) -M n j ) -ν d n (r+1)p-1 2 -d n rp-1 2 ,
and clearly

∆M n r - j∈Br M n j -u n+1 j ≤ 0, (31)
due to the correction step (21

)-(22). Then ∆M n r ≤ j∈Br M n j -u n+1 j that is 0 ≤ λ n r = ∆M n r j∈Br M n j -u n+1 j ≤ 1.
And finally

m n j ≤ u n+1 j = (1 -λ n r )u n+1 j + λ n r M n j ≤ M n j ,
International Journal on Finite Volumes since m n j ≤ u n+1 j ≤ M n j by construction. Thus (28) is proved for ∆M n r > 0. (The other case is proved using the same method.) It implies (29). It remains to prove (27). The definition of λ n r associated with equation (30) yields to

j∈Br u n+1 j = j∈Br u n+1 j + λ n r j∈Br M n j -u n+1 j = j∈Br u n+1 j + ∆M n r = j∈Br (u n j -ν(u n j -u n j-1 )) -ν(d n (r+1)p-1 2 -d n rp-1 2
).

Summing up with respect to r gives (27) and completes the proof.

Finite Volume Form

The next step is to identify the algorithm ( 2) to (26) as a standard finite volume scheme as stated in (35). For the simplicity of notations, assume that ∆M n r > 0. Due to the definition of the scheme (2) to (26), one has

u n+1 j = u n+1 j + λ n r (M n j -u n+1 j ) = ûn j -∆m n j + λ n r (M n j -u n+1 j ) = u n j -ν(u n j -u n j-1 ) -ν(d n j+ 1 2 -d n j-1 2 ) + -∆m n j + λ n r (M n j -u n+1 j ) .
In order to be able to rewrite the scheme as a finite volume scheme, we need to construct some fluxes e n j+ 1 2 such that for all j ∆t(e n j+ 1 2

-e n j-1 2 ) = -∆m n j + λ n r (M n j -u n+1 j ). ( 32
)
The solution can be constructed as

e n rp-1 2 = 0, ∀r ∈ Z, e n j+ 1 2 = 1 ∆t rp≤k≤j -∆m n k + λ n r (M n k -u n+1 k ) , rp < j. ( 33
)
This formula is correct because (24)-(26) implies that (33) is correct for j = (r + 1)p -1 : e n (r+1)p-1 2 = 0. If one (or more) of the ∆M n r < 0 is negative, the result remains the same.

Lemma 2.4

The L 1 norm of (e n j+ 1

2

) is bounded

||e n || 1 = ∆x j |e n j+ 1 2 | ≤ 2C 1 (2k + 1)p||u|| BV , (34)
where ||u|| BV is the BV norm of the solution.

Proof Due to the definition (33) we deduce

||e n || 1 ≤ p ∆t ∆x j |∆m n j |.
International Journal on Finite Volumes

But one always has

|∆m n j | ≤ |û n+1 j -(u j -ν(u n j -u n j-1 ))| ≤ ν(|d n j+ 1 2 | + |d n j-1 2 |) ≤ ν(|c n j+ 1 2 | + |c n j-1 2 |),
moreover due to the hypothesis (4) we finally get ∆x ∆t

j |∆m n j | ≤ 2C 1 (2k + 1) j |u n j -u n j-1 | ≤ 2C 1 (2k + 1)||u|| BV .
This ends the proof.

Then we are able to write the repair algorithm as

u n+1 j -u n j ∆t + a u n j -u n-1 j ∆t = (s n j -s n j-1 ), ( 35
)
where s n j is defined by

s n j = - 1 ∆x ad n j+ 1 2 + e n j+ 1 2 . ( 36
)
The L 1 norm of s n is bounded by

||s n || 1 ≤ C 1 (2k + 1)||u|| BV + ||e n || 1 ≤ C 1 (2k + 1)(2p + 1)||u|| BV . (37)
From (35) we get

u n+1 -u n ∆t + a ∆x (I -T )u n = (I -T )s n , ( 38
)
where T is the translation operator to the right. This is equivalent to say that

u n -(I + ν(T -I)) n u 0 = ∆t n-1 p=0 (I + ν(T -I)) n-p-1 (I -T )s p . (39)
The main convergence result is given by the Theorem 2.5 Assume that ν < 1. The difference between the upwind scheme and the repair algorithm (as defined by ( 2) to (26), or (35-36)) tends to zero with ∆x in the sense that there exists a constant C 2 > 0 which does not depend on k, p, C 1 , such that

||u n -(I + ν(T -I)) n u 0 || 1 ≤ C 2 C 1 (2k + 1)(2p + 1)||u|| BV (n∆t)∆x. ( 40
)
Proof Equation (38-39) shows that the scheme exactly fits in the framework developed in [START_REF]Després -Lax theorem and finite volume schemes[END_REF] and applied in [START_REF] Després | Proceedings of the Hispano-French Numerical Analysis School -Cadix[END_REF] for the convergence analysis of non linear schemes for linear advection. The key estimate proved in [START_REF]Després -Lax theorem and finite volume schemes[END_REF] is that

||(I + ν(T -I)) q+1 (I -T )|| 1 ≤ C 1 ν(1 -ν) 1 √ q . ( 41
)
International Journal on Finite Volumes

The constant C is a universal one. A simplified proof is provided in appendix for the case ν = 1 2 . This estimate means that the right hand side of (39) is made of terms that are small with respect to np. The final stage of the proof is straightforward. One has

||u n -(I + ν(T -I)) n u 0 || 1 ≤ ∆t max ||s p || 1 × n-1 p=0 ||(I + ν(T -I)) n-p-1 (I -T )|| 1 ≤ C 1 ν(1 -ν) ∆t max ||s p || 1 n q=1 1 √ q ≤ C 1 ν(1 -ν) ∆t max ||s p || 1 n 0 dy √ y ≤ C 1 ν(1 -ν) ∆tC 1 (2k + 1)(2p + 1)||u|| BV √ n 2 ≤ C 2a √ 1 -ν C 1 (2k + 1)(2p + 1)||u|| BV √ T ∆x. Defining C 2 = C 2a √
1-ν ends the proof. Remark: Since the numerical solution of the upwind scheme converges to the exact solution meaning

(I + ν(T -I)) n u 0 → u(n∆t) in L 1 (R) as ∆x → 0,
then the inequality of theorem 2.5 is a convergence result. Of course this inequality does not explain that the repair algorithm is better than the upwind scheme. Actually Repair is equivalent to the upwind scheme is one choose c n j+ 1 2 = 0, for all j and n. However the result shows that repair can not diverge if p 2 ∆x → 0.

Remark: On the other hand the error estimate blows up if p is too large. At least we need p 2 ∆x → 0 to get a vanishing error on the right hand side of the estimate. This estimate is the reason why we have incorporated the correction step in the repair algorithm.

Remark: One may wonder the reason of the 1 a √

1-ν contribution in the definition of the constant C 2 . Indeed if a → 0 or ν → 1, then C 2 can go to infinity which makes the estimate of convergence meaningless. First of all, for a given computation C 2 < ∞. Second of all, it is possible to use sharper estimates to get rid of the 1 a . This is done in [START_REF]Després -Lax theorem and finite volume schemes[END_REF] for the convergence analysis of TVD schemes for instance. Finally one may argue that a similar argument should be possible at the theoretical level to get rid of the 1 √ 1-ν . Nevertheless real computations use ν < 1 which is another reason to use this hypothesis. At the numerical level, we never saw any dependence of the rate of convergence with respect to this parameter. Thus we consider this as an artifact of the analysis.

Repair algorithm for remapping

In this section we show how to extend the previous result to take into account some very simple remapping algorithms. The main idea is that remapping is very close to transport. So let us consider a uniform mesh. (The uniformity of the mesh is not absolutely necessary however this hypothesis simplifies a lot the proof.). At time step n, the value of the unknown is u n j in the cell j. The boundaries of cell j are x n j-

t t n+1 n X j+1/2 X j+1/2 n X n j-1/2 X j-1/2 n+1 n+1 j u n j u n+1
1 2 and x n j+ 1 2 = x n j-1 2
+ ∆x, where ∆x is the mesh size. Then the mesh changes. We assume the simplest change: a uniform translation to the left, then

x n+1 j+ 1 2 = x n j+ 1 2 -a∆t, ∀j, (42)
where a > 0 is a kind of mesh velocity and ∆t is an equivalent time step. We assume that the mesh does not move too fast, that is we assume the CFL condition

ν = a ∆t ∆x ≤ 1.
A standard remapping algorithm is divided in two stages : first reconstruct a high order profile using a MUSCL like algorithm, then project on the new mesh. Let us detail these operations. An illustration is given in figure 1.

The reconstruction amounts to the definition of the slopes in function of the old values. So we define

v n j (x) = u n j + w n j (x), ( 43
)
where w n j is the slope function which is added to the average value at step n. A natural hypothesis is that the average value of w n j (x) is zero

x n j+ 1 2 x n j-1 2 w n j (x) = 0, (44)
and that w n j (x) is a local reconstruction, that is

|w n j (x)| ≤ C 3 j-k≤q≤j+k |u n q -u n q-1 |. ( 45
)
International Journal on Finite Volumes

Then we project the reconstructed function onto the new mesh. Thus

u n+1 j = x n+1 j+ 1 2 x n+1 j-1 2 v n (x)dx. ( 46
)
That is

u n+1 j = x n j-1 2 x n+1 j-1 2 v n j-1 (x)dx + x n+1 j+ 1 2 x n j-1 2 v n j (x)dx.
Replacing v n j (x) and v n j-1 (x) using (43), one gets

u n+1 j = νu n j-1 + (1 -ν)u n j + x n j-1 2 x n+1 j-1 2
w n j-1 (x)dx +

x n+1 j+ 1 2

x n j-1 2

w n j (x)dx.

Let us define

c n j-1 2 = - 1 x n j-1 2 -x n+1 j-1 2 x n j-1 2 x n+1 j-1 2 w n j-1 (x)dx = - 1 a∆t x n+1 j-1 2 +a∆t x n+1 j-1 2
w n j-1 (x)dx. ( 47)

Hypothesis (45) turns into

|c n j+ 1 2 | ≤ C 3 j-k≤q≤j+k |u n q -u n q-1 |. (48)
Moreover using (44) -the fact that the mean value of the correction is zeroyields to the definition of the new value u n+1 j

u n+1 j -u n j ∆t + a u n j + c n j+ 1 2 -u n j-1 -c n j-1 2 ∆x = 0. ( 49
)
With these notations the remapping algorithm is equivalent to a non linear transport. Thus we easily generalize the results of the previous section to remapping. The repair algorithm for remapping is then :

1) reconstruct a high order approximation as in (43-45) ; 2) make the correction step as in subsection 2.2 ; 3) repair.

The generalization of the correction step of subsection 2.2 is easy to write down, since the notations are quite similar between section 2 and section 3. It is sufficient to notice that once p ∈ N ⋆ , the size of the boxes, has been chosen, then the correction step checks if the inequalities (7-8) become true. If one of them is not true, then we need to minimize the value of c n j+ 1 2 at the borders of the boxes until the inequalities are true.

Then the convergence result of theorem 2.5 is extended to Repair remapping algorithms. For a given size of the boxes p, the difference between the Repair Algorithm and the upwind or donor cell method is bounded by estimate (40) (with C 1 = C 3 ).

International Journal on Finite Volumes

In the numerical result presented, we tried the Lax-Wendroff scheme plus Repair and the Downwind scheme plus Repair. The velocity is a = 1.

The Lax-Wendroff scheme is given by

(LW) u n+1 j -u n j ∆t + u n j + 1 2 (1 -ν)(u n j+1 -u n j ) -u n j-1 + 1 2 (1 -ν)(u n j -u n j-1 ) ∆x = 0.
The Downwind scheme is

(DW) u n+1 j -u n j ∆t + u n j+1 -u n j ∆x = 0.
The LW intends to be representative of high-order prediction schemes. The DW scheme intends to be representative of highly anti-dissipative prediction schemes.

Concerning the Repair algorithm we used different box sizes. We advected two initial profiles: a smooth one (u 0 (x) = -cos 2πx on [0, 1]) and a discontinuous one (u 0 (x) = 1 for 0 < x < 0.5 and u 0 (x) = 0 for 0.5 < x < 1) on a domain with periodic boundary conditions.

Figures

Essentially the results LW+Repair show good convergence properties independently to the size of the box p. See Figures 2 and3. However one sees for the cosine initial condition that p = 1 leads to some discrepancy at the extrema of the solution, that p = 50, 100 leads to a prediction of the extrema that is less accurate. On the other hand a relative small box of size p = 4 gives the best result. On this case, a small amount of local repairing gives good results. For the step initial condition of Figure 3, the spreading is quite the same with all box sizes. On the other hand the DW+repair results may be very sensitive to the size of the box. For a given number of cells, one gets the best results for the smaller size of the box, see Figure 4, staircases appear, exactly as with the UltraBee scheme as reported in [START_REF] Després | Lagoutière -Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF], see also [START_REF] Bouchut | An antidiffusive entropy scheme for monotone scalar conservation laws[END_REF] and [START_REF] Xu | Anti-diffusive flux corrections for high order finite difference WENO schemes[END_REF]. For a given box size the numerical solution converges to the exact one as the number of cells increases, see Figure 5. Unfortunately the global algorithm where the size of the box is equal to the number of cells (i.e. p = n) seems to diverge as the time increases, see Figure 6. Our experiments showed that this pathology is highly sensitive to the CFL number. For CF L = .5, the global algorithm seems to be correct. We retain that the global algorithm is not reliable for such test cases.

Order of convergence

In tables 1 to 3 are gathered some errors in various norms. The size of the boxes is p = 4. The final time is t = 1. We increase the number of cells from 100 to 1600. The order of convergence of the Lax-Wendroff plus Repair algorithm is ≈ 1 2 in L 1 for a discontinuous profile. This is in accordance with the theoretical result of Theorem 2.5. For a smooth profile the order is ≈ 2 in all norms.

International Journal on Finite Volumes Repair Algorithms in 1D ≈ 0 Table 2: Errors with respect to the size of the mesh. Lax-Wendroff plus Repair p = 4. The initial condition is u 0 (x) = 1 for 0 < x < 0.5 and u 0 (x) = 0 for 0.5 < x < 1. The final time is t = 1. The Courant number is CF L = 0.12345.

Compare with figure 3.

The order of convergence of the Downwind plus Repair algorithm is approximatively one in all norms for a smooth profile. It seems that Downwind plus Repair finds the closest step profile to the initial profile in a few time steps. Therefore Downwind plus Repair is almost exact for a step profile. This is very similar with what has been reported in [START_REF] Després | Lagoutière -Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF] for the Downwind TVD scheme (also known in the literature as the Ultra-Bee scheme).

For the discontinuous profile u 0 (x) = 1 for 0 < x < 0.5 and u 0 (x) = 0 for 0.5 < x < 1, the result are almost exact for 200 and more cells.

Gas dynamics in 1D

In this section we are solving the gas dynamics equations in 1D with a Lagrange+Remap code. This code is built on two components: a Lagrangian scheme and a Remap strategy which may need repair. The primitive variables are ρ the density, u the velocity and ε the specific internal energy, whereas, the conservative variables are m the mass, µ = mu the momentum and E = E + K = mε + 1 2 mu 2 the total energy (as the internal energy plus the kinetic energy). Density and specific internal energy are given on cell centers whereas velocity is given on nodes. The Lagrangian scheme is based on a compatible discretisation as described in [START_REF] Caramana | Shashkov -Elimination of Artificial Grid Distortion and Hourglass-Type Motions by Means of Lagrangian Subzonal Masses and Pressures[END_REF], [START_REF] Caramana | The Construction of Compatible Hydrodynamics Algorithms Utilizing Conservation of Total Energy[END_REF] and [START_REF] Campbell | Shashkov -A compatible lagrangian hydrodynamics algorithm for unstructured grids[END_REF]. This Lagrangian scheme is associated with a high-order remapping algorithm based on [START_REF] Loubère | Shashkov -A subcell remapping method on Staggered Polygonal Grids for Arbitrary-Lagrangian-Eulerian methods[END_REF]. This remapping algorithm is conservative in mass, momentum and total energy, reversible (if the old and new meshes are identical, the primitive variables remain unchanged) and high-order as a piecewise linear reconstruction is performed for each variable before remapping as [START_REF] Loubère | Shashkov -A subcell remapping method on Staggered Polygonal Grids for Arbitrary-Lagrangian-Eulerian methods[END_REF].

When solving the Euler equations we have to ensure that the density, the specific internal energy and the pressure remain positive during the calculation. Our Lagrangian scheme ensures this property whereas the remapping algorithm does not necessarily: a Repair method has then been added. Moreover the control of spurious oscillations is performed because a maximum principle is enforced by the repair method which was not an intrinsic property of the remapping method (see [START_REF] Loubère | Wendroff -The Repair paradigm: new algorithms and application to gas dynamics[END_REF]). The repair process has to be embedded into the remapping process as the primitive variables (ρ, u, ε) are not independent from the conservative variables (m, µ, E); for example momentum is computed using velocity (which may need repair) and mass. However the mass is computed using the density which may need repair as well. Therefore the Remap+Repair algorithm performs the following three steps:

1. Remap mass, then repair density -mass and density are final values; 2. Remap momentum, then repair velocity -momentum and velocity are final values;

3. Remap internal energy and kinetic energy, then repair specific internal energy -energies are final values.

The repair algorithm finally brings several new properties to the overall scheme, the respect of positivity and the preservation of a maximum principle.

It is well known that any numerical scheme which does not preserve the positivity by construction has to have a "kind" of non-conservative repair method such as ρ = max(0, ρ) and ε = max(0, ε) to exactly ensure positivity of density and specific internal energy. However we believe that a repair has to be performed more carefully especially to ensure the conservation. The problems chosen to test the method in 1D are the Sod Riemann problem, the Le Blanc Riemann problem and the Collela-Woodward blastwave. For each of these problems the code is run in Euler regime as Lagrange+Remap with a local repair method.

Sod Riemann problem

The Sod Riemann problem is defined on Ω = [0; 1] by two states of the same ideal gas γ = 1.4, separated at X = 1 2 : (ρ l , u l , p l ) = (1, 0, 1) and (ρ r , u r , p r ) = (0.125, 0, 0.1). The CFL is chosen equal to 0.25. In table 4 are gathered the convergence results in L p norms (p = 1, 2, ∞) for mesh sizes from 100 to 1600.

In Figure 8 are presented the density and the specific internal energy versus the exact solution for several mesh sizes 100, 200, 800. The Sod Riemann problem is a mild problem in the sense that the repair process can be canceled without generating any problem on the numerical solution. On the other hand this problem shows that a repair method, even not needed, is not polluting the numerical solution.

Le Blanc Riemann problem

The Le Blanc Riemann problem is defined on Ω = [0; 9] by two states of the same fluid γ = 5 3 , separated at X = 3: (ρ l , u l , p l ) = (1, 0, 2 3 10 -1 ) and (ρ r , u r , p r ) = (10 -3 , 0, 2 3 10 -9 ). The CFL is chosen equal to 0.25. In table 5 is gathered the convergence results in L p norms (p = 1, 2, ∞) for a mesh size from 180 to 2880. In this problem the specific internal energy simulation is the most difficult part because most of the numerical schemes (Eulerian or Lagrangian) converge very slowly on this problem. In Figure 9 are presented the energy for different mesh sizes (180,720 and 1440). The Le Blanc Riemann problem is a more severe test case where the repair method is necessary, otherwize negative internal energies are created and the code fails if no special treatment is provided. Moreover if only positivity is ensured then the overshoot observed in Figure 9 is much more pronounced.

Collela-Woodward blastwave

The Collela-Woodward blastwave problem simulates on Ω = [0, 1] the interaction of two Riemann problems defined by the three states (ρ l , u l , p l ) = (1, 0, 1000), (ρ m , u m , p m ) = (1, 0, 0.001), (ρ r , u r , p r ) = (1, 0, 100) of the same ideal fluid γ = 1.4 separated at X 1 = 0.1 and X 2 = 0.9. The CFL is chosen to be equal to 0.7. In figure 10 are presented the density for 200, 400, 1600 mesh sizes compared to a reference solution showing numerically the convergence of the method. This problem, as the Le blanc Riemann problem, can not run without a repair method. The repair process actually fixes the unphysical densities and specific internal energies and ensures

International Journal on Finite Volumes

Conclusion

In this paper we investigated the convergence in 1D of a repair algorithm. A Repair algorithm can be defined as a post-processing to conservatively ensure a maximum principle (hence the positivity) for any numerical scheme. In order to prove the main convergence result we introduced a box of size p in which the distribution of the excess of mass is performed. It seems that the size of the box (parameter p) is not necessarily an important parameter if one uses a high-order prediction scheme as the Lax-Wendroff scheme for non oscillating computations. But with a more antidissipative prediction scheme as the Downwind scheme, the results can vary with p. When p is too large the numerical solution may not be correct. It gives some indication that it is much preferable to restrict ourselves to local Repair (i.e with a local redistribution of the mass). An important feature of the Repair paradigm is the simplicity and versatility in any dimension. Moreover the repair process is independent of the kind of mesh used; cell-centered values or nodal values can be repaired the same way. We only need the notion of neighborhood to define the bounds and to redistribute the amount of conservative variable. Therefore any repair algorithm is suitable for staggered formulation where physical variables are not defined at the same place -classical Lagrangian schemes in CFD deal with cell-centered density, nodal velocity etc. The test cases in 1D show non only the convergence of the repair process on transport equation but on 1D gas dynamics as well. We have shown on several problems that the repair process is usefull to correct some bad behaviors of Lagrange+Remap schemes for a non-linear system of PDEs.

1 2≥ 0 and c n rp-1 2 < 0 : 2 < 0 and c n rp-1 2 ≥ 0 :

 120220 Then[START_REF] Loubère | Wendroff -The Repair paradigm: new algorithms and application to gas dynamics[END_REF] is true without condition. Thus ϕ 3,- Then it is not possible to simplify the inequality[START_REF] Loubère | Wendroff -The Repair paradigm: new algorithms and application to gas dynamics[END_REF]. Thus we impose ϕ 4,-

1 2< 0 and c n rp-1 2 ≥ 0 : 1 2.

 1201 [START_REF] Shashkov | Wendroff -The Repair Paradigm and Application to Conservation Laws[END_REF] is true without condition. Thus ψ 4,- each of the cases considered in inequalities[START_REF] Xu | Anti-diffusive flux corrections for high order finite difference WENO schemes[END_REF] to (17). We gather the restrictions it imposes for all µ n rp-The mathematical definition of the correction algorithm is the following. Definition 2.1 Let us define the corrected fluxes at the boundaries of the boxes

Figure 1 :

 1 Figure 1: Remapping stage. The old and new densities are the bullets. The high order reconstruction stage are represented with the slopes.

Figure 3 : 1 International

 31 Figure 3: Lax-Wendroff plus Repair : 100 cells. The initial condition is u 0 (x) = 1 for 0 < x < 0.5 and u 0 (x) = 0 for 0.5 < x < 1. The size of the box is p = 1, 2, 4, 10, 50, 100 from top left to bottom right. The final time is t = 1

Figure 4 : 1 InternationalFigure 5 : 1 International

 4151 Figure 4: Downwind plus Repair. The initial condition is u 0 (x) = -cos 2πx. The size of the box is p = 1, 2, 4, 10, 50, 100 from top left to bottom right. The number of cells is N = 100. The Courant number is CF L = 0.12345. The final time is t = 1

Figure 6 :

 6 Figure6: Downwind plus Repair : 100 cells. The initial condition is u 0 (x) = 1 for 0 < x < 0.5 and u 0 (x) = 0 for 0.5 < x < 1. The size of the box is p = 1, 2, 4, 10, 50, 100 from top left to bottom right. The final time is t = 1. The results are quite good, except for the global algorithm which seems to be in advance.

Figure 7 :

 7 Figure 7: Downwind plus Repair : 100 cells. The initial condition is u 0 (x) = 1 for 0 < x < 0.5 and u 0 (x) = 0 for 0.5 < x < 1. The size of the box is p = 100. The final time is t = 1, 2, 3, 4, 5, 6. The global algorithm is in advance.

Figure 8 :

 8 Figure 8: Sod Riemann problem -Euler regime (as Lagrange+Remap) with Repair process -Density for a 100, 200, 800 mesh sizes vs exact solution.

Figure 9 :

 9 Figure 9: Le Blanc Riemann problem -Euler regime (as Lagrange+Remap) with Repair process -Specific internal energy for a 180, 720, 1440 mesh sizes vs exact solution -Full size (left) and zoom (right).

Figure 10 :

 10 Figure 10: Collela-Woordward blastwave problem -Euler regime (as La-grange+Remap) with Repair process -Density for 200, 400, 1600 mesh sizes vs reference solution -Full size (left) and zoom (right).

Table 1 :

 1 Errors with respect to the size of the mesh. Lax-Wendroff plus Repair p = 4. The initial condition is u 0 (x) = -cos 2πx. The order is approximatively 2. The Courant number is CF L = 0.25. Compare with figure 2.

	International Journal on Finite Volumes

Table 3 :

 3 Errors with respect to the size of the mesh. Downwind plus Repair. The initial condition is u 0 (x) = -cos 2πx. The order is approximatively 1. Compare with figure5.

	cells	L 1	L 2	L ∞
	100	0.090947985213994 0.115785613616216	0.2710309
	200	0.041511107024698 0.051893396383501 0.124615354
	400	0.023674021241459 0.031497599561389 0.101446636
	800	0.012277712311722 0.016193499317599 0.052800138
	1600 0.006085189480290 0.008036148646353 0.029193723
	order	≈ 1	≈ 1	≈ 1

Table 4 :

 4 Errors for density with respect to the size of the mesh. Sod Riemann problem -Euler regime (as Lagrange+Remap) with Repair process.

	International Journal on Finite Volumes

Table 5 :

 5 Errors for specific internal energy with respect to the size of the mesh -Le Blanc Riemann problem -Euler regime (as Lagrange+Remap) with Repair process.

	International Journal on Finite Volumes

International Journal on Finite Volumes

For gas dynamics, the maximum principle means that if m, µ, E are remapped and ρ, u, ε are the associated repaired variables, then no new ewtrema is created in the solution ρ, u, ε of this system. For a general system it is an open problem to design a criterion for repairing International Journal on Finite Volumes

Acknowledgments

The authors thank M. Shashkov and B.Wendroff for fruitful discussions. This work LAUR-04-7960 was performed as the first ans second authors were respectively an invited visitor and an employee of the Los Alamos National Laboratory during summer 2004. This work was performed under the auspices of the US Department of Energy at Los Alamos National Laboratory, under contract W-7405-ENG-36. The authors acknowledge the partial support of the DOE/ASCR Program in the Applied Mathematical Sciences and the Laboratory Directed Research and Development program (LDRD). The authors also acknowledge the partial support of DOE's Advanced Simulation and Computing (ASC) program.

A A simple proof of inequality 41 in the case ν = 1 2 We prove here that there exists C > 0 such that for ν = 1 2 one has

A complete proof is in given [START_REF]Després -Lax theorem and finite volume schemes[END_REF], that also takes into account variable time steps. One has

The binomial coefficient is q j = q! j!(q-j)! . The series j → q j is increasing from 0 to [ q 2 ] ± 1 and is decreasing from [ q 2 ] ± 1 to m, where [ q 2 ] ∈ N is the entire part of q 2 . Thus

For a positive number m, the Stirling formula [START_REF] Abramowitz | Handbook of mathematical functions[END_REF] gives that m! ≈ √ 2πm( m e ) m . Plugging in the right hand side of our last inequality, one gets that

From (51-53) one gets

It ends the proof of the lemma in the case ν = 1 2 . This proof can easily be generalized for ν = 1 2 .

International Journal on Finite Volumes


Claude Vanhecke

consultant. L A Assouère

P Wang

F Durand-Estèbe

J M Caignet

M Dilhac

Bafleur

Anqing Wang

Laurent Assouère

Paul Durand-Estèbe

Fabrice Caignet

Senior Member, IEEE Jean-Marie Dilhac

THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <

Keywords: Energy harvesting, energy scavenging, structural health monitoring, ultra low-power converter, multisource harvesting, battery-free, nanowatt voltage reference

We suggest an innovative architecture for an efficient energy generator devoted to the powering of a wireless sensor network deployed for Aircraft Health Monitoring. This battery-free generator captures energy from its environment (transient thermal gradients as a main source, and vibrations as a secondary source allowing early biasing of the generator) and stores this energy in ultracapacitors. In this way, this multi-source architecture benefits from the synergy between energy scavenging and harvesting: vibrations bring low but early and permanent energy. They also contribute to energy harvesting during cruise while thermal gradients have vanished. The use of active diodes and of a very low bias current of 10nA/branch allow achieving ultra-low power consumption, experimentally demonstrated on two different CMOS technologies. It is also proven that enough energy could be delivered to power the functions of a wireless sensor node.

and of impact damages taking place during the aircraft life. It is generally considered as a powerful tool to decrease inspection costs, to optimize margins in mechanical design, and consequently, to reduce aircraft weight, fuel consumption and emissions of greenhouse gases.

Only limited implementations of AHM have already been done using wired technology such as Acoustic Airframe Monitoring

System by Ultra Electronics Ltd [START_REF]Ultra Electronics Aircraft systems[END_REF]. However, because of the gain associated with a wireless architecture, support of AHM by Wireless Sensors Networks (WSN) based on Micro Electrical Mechanical Systems (MEMS) technology, is envisioned by major aircraft manufacturers.

In the above context, a network of a relatively large number of self-powered MEMS-based nodes would perform sensing, data processing and wireless transmission of information. However, the main technological barrier to the wide development of WSN is their energy autonomy.

In such a WSN, the most obvious wireless energy supply system consists of primary batteries. However, whatever the batteries performance and size, they only store a limited amount of energy and exhibit a limited lifetime therefore placing an unacceptable upper limit on the network lifetime itself (given the fact that nodes may be placed in remote areas with very limited access). As a result, the replacement of hundreds of dead batteries would induce a prohibitive maintenance cost together with an environmental issue for their disposal. More important, at temperature levels encountered at high altitude in unpressurized areas far from engines (typically -60°C) their efficiency is drastically reduced whereas unacceptable safety issues (thermal runaway and fire) are raised. Fortunately, primary batteries can be eliminated through the use of environmental energy capture techniques [START_REF] Mitchell | Energy Harvesting Applications and Architectures at Boeing Commercial Airplanes[END_REF], which use an energy conversion transducer tied to an integrated rechargeable power storage device, then enabling the wireless sensor node an almost infinite lifetime.

In this paper, we present an ultra-low power converter for a multisource battery-free energy generator dedicated to aeronautics applications that would enable almost infinite energy-autonomy to a WSN node. It provides a regulated voltage with a very low quiescent current consumption (≤300nA). The proposed architecture is based on two energy sources and on ultracapacitors for storage. The paper is organized as follows:



Section II presents the potential energy sources available for the targeted SHM application and in particular, the two selected sources: thermal gradient and mechanical vibrations.

 Section III discusses the different design choices for the architecture of the multisource energy generator.

 Section IV provides the experimental results obtained on two different CMOS technologies and some comparison with commercially available solutions.

II. ENERGY GENERATION IN AERONAUTICS ENVIRONMENT

For energy capture, two principles may be considered, called energy harvesting (continuous source) and energy scavenging (intermittent source). However, availability of energy is in both cases limited, and for such a self-powered network, energy is therefore a critical issue, and hardware design must consider energy as a main constraint.

It is worth mentioning that coupled energy sources can be considered since it may occur that a single category would not be enough to power a node during all phases of a flight. More precisely, time shifts in the availability of environmental energy together with intrinsic different time constants of the transducers and the by-nature synergy between scavenging and harvesting may praise for such a multi-source configuration.

Various types of ambient energy may be considered:

-photovoltaic, in the very limited case of outside sensors and daylight flights [START_REF] Meekhun | Design of a solar harvester system for a wireless sensor network deployed for large aircraft in-flight tests[END_REF] -wind or acoustic energy, which we are considering in another on-going work dealing with aeroacoustic [START_REF] Montheard | Proof of concept of energy harvesting from aero acoustic noise[END_REF], with the drawback of the required development of very specific harvesters -low (compared to acoustic) frequency vibrations -thermal gradients

The two last energies are of interest provided that the operating principle and packaging of the harvesters are compatible with the harsh environment, mainly in terms of temperature and acoustic pressure when implemented in the engine area. Worth to mention is also the strong governing principle of aeronautics: volume and weight are the ennemies. In other words, if the harvesting system is much heavier than the cable that would power the sensor in a wired configuration, the chance of accepting ambient energy harvesting is low.

Although in terms of energy loss, immediate use of the captured environmental energy by the WSN node would be more efficient, an energy-storage device is required as an energy buffer between the WSN node and the energy source. Additionally, long-term energy storage may be desired to budget for future energy consumption when scavenging or harvesting efficiency is low.

In the aeronautics application context, the use of secondary (rechargeable) batteries is prohibited, as they suffer from even worse environmental limitations than primary ones [START_REF] Meekhun | Charge and Discharge performance of Secondary Batteries according to Extreme Environment Temperatures[END_REF]. Electrical double-layer capacitors, also named super-capacitors or ultracapacitors, are the solution for transient storage. They store electrostatic energy between a solid electrode and oppositely charged electrolyte ions. They offer a high capacitance in a small volume together with a virtually infinite lifetime. However, conversely to batteries, they do not provide a fairly constant output voltage, the output of an ultracapacitor dropping from full value to zero during discharge: a voltage regulator is therefore required for biasing the signal processing part of the WSN node. The ultracapacitor from Maxwell we have selected had been extensively tested vs. temperature in a previous work. We had experimentally demonstrated that its performance was unaffected down to -55°C [START_REF] Meekhun | Design of a solar harvester system for a wireless sensor network deployed for large aircraft in-flight tests[END_REF].

It must be stressed that, at this stage, the parameter of importance is the transfer through the transducer of a maximum of energy from the environment to the ultracapacitor, whatever the conversion efficiency, all energy not captured being ultimately lost. Conversely, energy treatment such as voltage regulation must favour efficiency.

In the following, the considered AHM applications preclude external devices mounted on the airfoils outside the aircraft, therefore excluding the use of solar energy, which is known to be one of the most efficient (in terms of power density) way of capturing energy from a system environment. Similarly, for deployment we consider non-pressurized areas, far from any heatgenerating equipment.

Since from take-off to cruise altitude, atmospheric temperature usually strongly decreases, thermal gradients may therefore appear inside and around the WSN node hardware. Referring to the International Standard Atmosphere (ISA) defined by the International Civil Aviation Organization (ICAO), at an en-route flight altitude of 12000 m, atmospheric temperature is about -60°C [START_REF]Standard Atmosphere[END_REF]. The above gradients can be converted into voltage via a thermoelectric generator (TEG) based on the Seebeck effect. However, sooner or later, the node and its surroundings will be back to an isothermal state and no more thermal flux will take place. Therefore, to both increase the gradient value and duration, we chose to capture the transient thermal gradient between aircraft mechanical structure and a phase change material [START_REF] Bailly | Energy Scavenging based on Transient Thermal Gradients: Application to Structural Health Monitoring of Aircrafts[END_REF][START_REF] Samson | Energy Harvesting for Autonomous Wireless Sensor Nodes in Aircraft[END_REF][START_REF] Kiziroglou | Heat storage power supply for wireless aircraft sensors[END_REF]. We identified water as a pertinent choice, therefore benefiting from its high specific heat (4.2 J/g.°C -to be compared with that of Si, 0.7 J/g.°C). Water freezing point (0°C) being situated within the aircraft temperature usual operating range, the system also benefits from energy exchanges related to water latent heat of fusion (330 J/g), with the added advantage of water being maintained at a constant temperature during the water freezing or ice melting phases, hence increasing temperature gradient values. Fig. 1 shows a cross-section drawing and a photograph of the experimental setup that allowed validating the thermal-gradient energy-scavenging principle. To roughly assess the upper limit for the thermal energy that could be captured, let us consider a sealed tank containing 1 g of water. As already mentioned, water exhibits a specific heat C of around 4.2 J/g.°C when liquid, 2.1 J/g.°C when solid, and an enthalpy of fusion H f of 330 J/g. Stating a water capsule temperature at take-off of 15°C, and a cruise temperature of -60°C, the energy W T associated with temperature variation ΔT and phase change is given by:

W T = ΔT solid .C solid + ΔT liq .C liq + H f = 519 J/g (1)
This value may be doubled as temperature variations take place both during climb and descent that is a total of 1038 J/g. This value may be slightly higher if take-off temperature and therefore initial temperature of the sensor node is higher. Unfortunately, it can conversely be very limited if this temperature is below icing temperature. As shown in Fig. 2, thermal energy can be converted into electrical energy by thermoelectric generators using the Seebeck effect. Such generators are current sources converting temperature difference between a cold and a hot surface into electrical energy. Despite the efficiency of such components that may be as low as 1% or less, the energy available at the TEG output is around 15J. This efficiency could be even more improved using an appropriate phase change material and an optimised heat storage unit [START_REF] Kiziroglou | Heat storage power supply for wireless aircraft sensors[END_REF].

Although the above method is attractive, energy scavenging will only start once the aircraft is climbing; this may raise operational issues if energy-storage devices are initially empty. This latter situation may often occur since ultracapacitors are known to suffer from large self-discharge currents. In particular, monitoring of the structure during take-off is of great interest since it is submitted to high strains. To cope with this drawback, we devised a system with a secondary energy source harvesting the mechanical energy associated with the vibrations of the aircraft structure. In an airliner, vibrations are maximum in engine area, and conversely are much less in cabin, with the exception of low (less than 10 Hz) accelerations associated with gust or turbulence. Fig. 3 qualitatively shows the experimental vibration spectrum recorded in close proximity of the jet engine of a large airliner. It is a multi-band spectrum, made of peaks regularly spaced. The fundamental frequency corresponds to the rotating speed of the fan of the engine low-pressure compressor. Its value -and those of the harmonics -are consequently depending upon flight phase, i.e take-off, climb, cruise... Moreover, the frequency peak is only a few hertz wide. These two characteristics make energy harvesting a difficult task when using narrowband piezo harvesters. To harvest vibrations, we devised a MEMS harvester based on a seismic mass applying the vibrations to piezoelectric beams [START_REF] Durou | Micromachined bulk PZT piezoelectric vibration harvester to improve effectiveness over low amplitude and low frequency vibrations[END_REF] [START_REF] Durou | Vers l'autonomie énergétique des réseaux de capteurs embarqués : conception et intégration d'un générateur piézoélectrique et d'un micro dispositif de stockage capacitif en technologie silicium[END_REF]. It delivers an alternating output which spectrum exhibits a fundamental frequency around 60Hz, i.e. at the rotating speed of the engine inlet fan. These vibrations originating from aircraft engines are not intense enough to power all the functions of the WSN node by themselves. However, as will be shown later, they are sufficient for biasing the electronics associated with the TEG alone as soon as the engines start, therefore saving time with respect to the scavenging of thermal gradients alone.

III. ARCHITECTURE OF THE MULTISOURCE BATTERY-FREE ENERGY GENERATOR

As already mentioned, the large deployment of wireless sensors, in particular for SHM is still limited by the issue of energy autonomy of the sensor node. Energy harvesting is an attractive solution to provide battery lifetime extension or even almost infinite autonomy to a battery-free sensor node. In this context, many authors have contemplated multisource harvesting to achieve a generic power generator for wireless sensors [13] [14]. In [START_REF] Colomer-Farrarons | A Multiharvested Self-Powered System in a Low-Voltage Low-Power Technology[END_REF], four sources (solar, vibrations, RF, thermal gradients) are combined and in [START_REF] Lhermet | Efficient Power Management Circuit: From Thermal Energy Harvesting to Above-IC Microbattery Energy Storage[END_REF], thermal energy harvesting is coupled to RF power to charge a microbattery. Both papers mainly focus on the validation of multiharvesting concept and its related power management. In both cases, the storage unit is placed after the voltage regulation block. Even more complex multi-source energy platforms have been theoretically envisioned in the literature accelera&on) frequency) fundamental*frequency* [START_REF] Beigné | Bringing robustness and power efficiency to autonomous energy harvesting microsystems[END_REF]. They not only deal with energy harvesting, but perform sensing and communication management from the point of view of energy transfer and storage efficiency. Despite the fact that these functions are ultimately mandatory in a wireless sensor, we have not implemented them in the more limited context of our work.

The main originality of the architecture we propose for a multisource energy generator (Fig. 4), consists in the combination of complementary energy sources for a rapid start-up of the self-powered battery-free system. It was initially motivated by aeronautics applications but is also applicable to other SHM scenarios and is compatible with other energy sources. The basic principle is to use a TEG as the primary source of energy whereas a piezoelectric generator is used as a secondary source. The latter one provides the necessary bias for the rectifier circuit required for the TEG and partially compensates the WSN node energy consumption during cruise over a long flight. The energy provided by the secondary source, through a voltage doubler with external Cin capacitor, is stored after rectification in a small capacitor Cbias that biases the voltage reference circuit (PTAT/BGR block). This circuit provides the bias current for the amplifiers in the active diodes. Once activated, these diodes allow harvesting energy from the TEG and storing it in ultracapacitors (UC1 and UC2) to be chosen with optimized leakage currents to minimize losses. To simplify circuit design, we have not considered implementing a Maximum Power Point Tracking (MPPT) procedure such as the architecture used in [START_REF] Lefeuvre | Buck-Boost Converter for Sensorless Power Optimization of Piezoelectric Energy Harvester[END_REF]. Consequently, depending upon ultracapacitors charge level, power captured is not permanently maximized.

Finally, voltage regulation is implemented using a low-dropout (LDO) regulator based on a PMOS power transistor. The energy generator should be able to supply a wireless sensor node. To allow an adequate use of the available energy, a careful power management of the WSN node (appropriate periodic measurement and communication, sleeping mode management) has to be implemented to save energy.

The main design challenges for the energy generator are:



Firstly, maximizing the transfer of energy from the transducer to the storage devices.



Secondly, as the voltage on the ultracapacitor storage devices will vary according to their discharge, a voltage regulation is needed. In this case, priority should be given to conversion efficiency.

 For both of these design challenges, another important requirement is microscale compatibility thus prohibiting the use of bulky passive devices as the ones needed for some impedance matching strategies [START_REF] Gyorgy | Review of Power Conditioning for Kinetic Energy Harvesting Systems[END_REF], for example.

 Finally, harvested energy being limited, the quiescent power consumption of the energy generator should be as low as possible. This is required for two reasons: the first one is that many SHM scenarios are using duty cycles with much longer periods in quiescent mode than in the active one and the second is related to the self-discharge current of ultracapacitors which is in the order of µA. The quiescent current of the power converter should not be higher than this value or even be made negligible compared to it. A trade-off will then have to be made between efficiency and power consumption.

For this latter challenge, we decided for the design of each analog block of the circuit to favour solutions minimizing the number of transistors.

A. Energy transfer maximization

As the energy harvested through thermal gradient results in two voltage polarities, even if there is only one polarity switching, a rectifier circuit is mandatory before storage. The classical solution consists in using a diode bridge such as the Graetz circuit, Fig. 5(a). The first drawback of this solution is that it does not allow maximizing the transfer of energy. For example, if the first voltage alternation is higher than the second one (it is the case in Fig. 2), it will not be possible to harvest the energy from this second alternation. To cope with this problem, the solution consists in implementing two ultracapacitors for the storage as shown in Fig. 5(b): UC1 for the positive alternation and UC2 for the negative one.

The second drawback of the Graetz circuit is that two diode thresholds are lost for the rectification of each voltage alternation.

To avoid this power loss, we then implemented the Graetz circuit using active diodes [START_REF] Gyorgy | Review of Power Conditioning for Kinetic Energy Harvesting Systems[END_REF][START_REF] Dallago | Comparison of Two Autonomous AC-DC Converters for Piezoelectric Energy Scavenging Systems[END_REF][START_REF] Sun | An Integrated High-Performance Active Rectifier for Piezoelectric Vibration Energy Harvesting Systems[END_REF] as shown in Fig. 5(c). However, as shown in Fig. 6, the active diode requires a voltage reference circuit to bias the amplifier (g n and g p bias signals provided by PTAT block in Fig. 5(c)). The design of this amplifier is critical and three important properties are required:  Ultra-low power consumption to comply with our requirement of a low quiescent current of the full converter.



Threshold voltage of MOS transistors as low as possible and independent of process variations as well as layout mismatching. The main issue could be the oscillation of the driving signal provided by the operational amplifier to the gate of the MOS switch. For example, in the 0.35µm HV CMOS technology, we used isolated high-voltage (20V) MOS transistors with a typical threshold voltage of 0.5V.

 High DC gain.

This latter condition, although requiring large size transistors, guarantees the second condition.

For the start-up of the active diodes, a minimum voltage around 0.9V (V th + V Dsat ) is required to bias the current mirrors and the differential pairs shown in the electrical schematics of Fig. 6. Once this voltage is reached, the threshold voltage of the active diode is given by the source to drain voltage of the MOSFET switch. Thanks to the regulation loop introduced by the operational amplifier that is designed with a high DC gain, this voltage is equal to the offset of the operational amplifier and can then be very small (few tenths of mV are targeted) compared to a classical diode threshold voltage. Nevertheless, before reaching this start-up voltage, the parasitic body diode of the MOSFET switch allows starting the rectifier with a standard diode threshold (~ 0.7V). In the case of totally discharged ultracapacitors and a single energy source (namely TEG), several minutes are needed to reach 0.9V.

B. Energy efficiency

We optimized the design of the active diodes with the objectives to get a threshold voltage in the range of 10 mV whereas minimizing the quiescent current consumption in the range of few hundreds of nA. To generate such a low bias current (10nA/branch), we designed a nanowatt voltage and current reference [START_REF] Vanhecke | Convertisseur à faible consommation pour la récupération d'énergie ambiante combinant deux sources pour application aéronautique[END_REF][START_REF] Rossi | Ultra-low Power CMOS Cells for Temperature Sensors[END_REF][START_REF] Vanhecke | Current generator, notably for current of the order of nano-amperes, and voltage regulator using such a generator[END_REF]. To this purpose, we combined subthreshold MOS devices with transistors working in the linear region and did not use any resistor as shown in Fig. 7. In this circuit, all transistors are working in saturation except N1 and N2 NMOS transistors that are in weak inversion (subthreshold region) and N3R, a very long channel transistor operating in the linear region, which acts as a high-value resistor [START_REF] Vanhecke | Current generator, notably for current of the order of nano-amperes, and voltage regulator using such a generator[END_REF]. Since N1 and N2 are working in the subthreshold region, the source voltage of N1 and then the drain voltage of N3R are given by:

!"! = !"!! = ! ln ! !! ! !! ! !! ! !!
where is ! = is the thermal voltage, proportional to absolute temperature (PTAT), and S N1 , S N2 , S P1 , S P2 are W/L ratios of the MOSFET's.

As a result, the current flowing in N3R transistor is:

= ! ! ! !"#!! ln ! !! ! !! ! !! ! !! (3)
with R onN3R , the on-resistance of N3R transistor. This current is then proportional to the absolute temperature. To define a very low reference current (~10nA), a high value resistor (~5MΩ) is needed. To reach such a value, the channel length of N3R

transistor was designed as long as 330µm. It has to be noticed that, for matching reasons, N4 transistor has also the same very long channel.

The expression of the reference voltage V ref , that is equal to the gate voltage of N4, is the following:

!"# = 2 !! + !! (4)
where !! is the transconductance of transistor N4 and V th its threshold voltage. Since V th has a negative temperature coefficient, V ref can be made independent of temperature. The value of V ref can then be adjusted through the size of the transistors. We have chosen to define it at 0.8V. In the 0.8µm SOI CMOS technology, respective aspect ratios W/L in µm of the transistors are: 60/30 for P1, P2 and P3, 20/20 for N1 and 160/20 for N2, 6.4/330 for N3R and 3.2/330 for N4.

A very similar current source subcircuit was proposed by Ueno et al [START_REF] Ueno | A 300 nW, 15 ppm/ C, 20 ppm/V CMOS Voltage Reference Circuit Consisting of Subthreshold MOSFETs[END_REF] to generate a nanowatt voltage reference. The main differences are the choice of using all the transistors in the subthreshold region except for the long-channel transistor, a more complex bias voltage subcircuit and the need for a start-up circuit. Compared to our proposed voltage reference, it has better performance regarding temperature sensitivity (15 ppm/°C vs 290 ppm/°C) and line sensitivity (20 ppm/V vs 800 ppm/V).

Nevertheless, the reduced number of transistors of our proposed voltage reference allows a very low quiescent power consumption of 600 nW with a 3V power supply for the whole multisource converter to be compared to the 300 nW at 1.5V for the single voltage reference circuit. Cbias capacitor when using only a TEG for energy generation ("Vbias without piezo") and the combination of a TEG with a piezoelectric transducer ("Vbias with piezo"). Also plotted is the TEG output voltage, Vs TEG, used for the simulation. Voltage output of the piezoelectric transducer, Vpz, is an AC voltage with ±3V amplitude and 60Hz frequency.

In addition, to limit the effects of process mismatches, particular care was taken for the layout and routing of the amplifier and voltage reference circuit (use of dummy devices, common centroid patterns...). Moreover, the very long channel transistors, N3R and N4, are designed in an isolated well including guard rings and good well contacts to insure a constant biasing of the substrate along the long channel.

As mentioned earlier, the originality of the proposed energy generator circuit is the combination of the TEG with a piezoelectric transducer to provide a rapid start-up of the energy scavenging from the TEG and a complementary source of energy during cruise. This beneficial effect is demonstrated on the simulation results of Fig. 8 that compares the case of the energy generation using the TEG alone with the TEG combined with the piezoelectric transducer. It can be noticed that using a multisource approach, the energy generator is operational as soon as the engines start.

Another important parameter that allows maximizing the energy transfer is the value of the storage ultracapacitor (UC). The time required to charge an UC depends on the time constant R S .UC of the circuit. This is particularly relevant since the energy generation is only transient during take-off and descent. If the value of the UC is too small, it will be rapidly charged but its maximum operating voltage, together with its small value, will limit the stored energy. On the contrary, if its value is too large, the time constant will prevent collecting the maximum of energy since the UC will not have the time to reach the voltage saturation. To optimize this trade-off, we simulated the charge of an UC via an ideal diode during the aircraft take-off and compared the open-circuit voltage to the one of a matched load. In our case, given the fact we use two UCs in series, the optimized value of each one is 1F [START_REF] Vanhecke | Convertisseur à faible consommation pour la récupération d'énergie ambiante combinant deux sources pour application aéronautique[END_REF]. It has to be noticed here that the piezoelectric transducer can also directly provide energy to the output, the active diode Dp3 providing isolation from the UCs since the piezoelectric transducer would not be able to directly charge them given the small current (~µA) generated by mechanical vibrations in the aeronautics environment. However, this direct connection to the output V OUT could be very useful to provide energy to the WSN node either on ground, as soon as aircraft engines start, or during long flights (several hours) since in these cases, the TEG might be inactive. Finally, a voltage regulation stage is required after the UC. Ideally, it should exhibit high conversion efficiency. The best candidate would be a switching regulator. However, given its complexity, in terms of quiescent power consumption, it will be less efficient than a simple linear regulator. We then chose to implement a low-dropout regulator (LDO) using a PMOS power transistor because of its simpler driving (no charge pump needed). The proposed nanowatt bandgap reference (BGR) circuit based on MOS devices allows providing the voltage reference for the LDO that is designed to supply a maximum peak current of 20mA.

The proposed energy generator provides two regulated voltage values: Vref or 2Vref that can be switched using an external command (Fig. 9). It has to be mentioned that in this circuit, N3R, N4 and N6 transistors have both a very long-channel and can be designed in the same isolated well.

Regulated voltages are required for supplying signal processing and memory circuits. For less sensitive circuits, a direct power supply from the UCs could be implemented. This would have the advantage to greatly increase the global efficiency.

IV. EXPERIMENTAL RESULTS

This multisource and battery-free energy harvesting architecture was validated on two technologies: a high-voltage 0.35µm CMOS technology from AMS available via the French Multi-Project Chip (CMP) service and a smart power 0.8µm Bipolar CMOS DMOS merged technology on SOI (TFSMART1) provided by Telefunken Semiconductors [START_REF]SMARTIS 1, Design Manuel[END_REF]. Optical views of these test chips are given in Fig. 10. On the first test chip (0.35µm CMOS), we validated the nanowatt voltage reference and the global architecture with a single source, i.e. the thermogenerator. On the second one (0.8µm SOI), we validated the multisource architecture. Both technologies provide high-voltage devices. A high-voltage technology is required since both the chosen thermoelectric generator, a Micropelt GmBH (MPG-D751), and the laboratory-made piezoelectric generator [START_REF] Durou | Micromachined bulk PZT piezoelectric vibration harvester to improve effectiveness over low amplitude and low frequency vibrations[END_REF] [12] could generate voltages as high as 10V. Both technologies allow isolating the active diodes and the very long channel transistor. The SOI technology has the additional advantage to completely eliminate any parasitic structure that could be triggered on in severe environment.

A. 0.35µm CMOS testchip

In this technology, several types of transistors are available. For the needs of this circuit, we chose transistors with the lowest threshold voltage and with an annular layout to limit leakage current and improve radiation immunity. Two versions of the circuit were designed on the same die: an ultra low-power (ULP) with 10nA/branch and a low-power (LP) with 100nA/branch.

The total size of the chip is 1.7mmx1.7mm.

Both LP and ULP versions of the circuit are functional. We only report here the results for the ULP circuit. Fig. 11 describes the block diagram of the ULP circuit. In addition to the functional pins, a dedicated pin called Temp allows sensing the silicon die temperature via a diode. Vbgr pin allows checking the functionality of the bandgap reference circuit as well as implementing an output voltage of 0.8V by short-circuiting it to Vout pin. Vin+ and Vin-are the pins for the connection of the thermogenerator (TEG). Vpol is the pin for the connection of the small capacitor Cbias that has a value of 10µF. This capacitor allows the selfbiasing of the circuit. Three pins are implemented for the ground rails: two for the internal circuit and one for power (GndP).

Vst+ and Vst-are the pins for the connections of the two ultracapacitors, UC1 and UC2 whose value is 0.9F. As already mentioned, the circuit allows providing two values for the regulated voltage Vout, which can be chosen by externally connecting the pins Vout and Vbgr. The targeted values were 1.6V and 0.8V, respectively. This could also be implemented via a specific internal circuitry triggered by an external logic signal.

For the targeted application of AHM that requires a high-energy autonomy, it is very important to choose an ultracapacitor with the best long-term self-discharge properties. These properties depend on the way the ultracapacitor was previously charged.

Among commercially available devices, Maxwell HC series, after being charged with a constant 5mA current, presents a selfdischarge of only a few % after 12h that is fully compatible with our application. This self-discharge is even smaller if the charge current is lower than this value, which will be mainly the case for energy harvesting applications. In the same way, the storage capacitor Cbias (10µF) also needs the same property of low leakage current to guarantee a good efficiency of the energy generator.

For the functionality testing, the energy generator is emulated using an electrical-waveform generator. First of all, we checked the functionality of the active diodes. As shown in Fig. 12, both diode types provide a very low diode threshold around 20mV, as expected from the simulations.

We also measured the values of the regulated voltage and of the quiescent current and their dispersion over the 10 dies received from the CMP service. As shown in Fig. 13, the measured values perfectly fall within the range defined by Monte-Carlo simulations. We also checked the temperature behaviour of the circuit. The temperature specifications of the application covers a wide range from -50°C to 100°C. As can be seen from Fig. 13, the temperature performance is a variation of 290 ppm/°C that is sufficient for the application and quite good for a PTAT/BGR block only consuming 100nA and without any offset compensation circuitry. For the validation of the energy generation from transient gradients, we were not able to perform in-flight testing. However, our partner, Airbus, provided typical temperature profiles that we used in our experimental test bench.

To simulate the aircraft temperature conditions, we tested the TEG scavenging module of Fig. 1 into an environmental test system TP04200A ThermoStream® from Temptronic corporation. To model the evolution of the temperature during a flight, we chose the standard atmosphere model of ICAO, that is a temperature of 15°C (59°F) and a pressure of 1013.25 hPa (14.7 psi) at sea level, and a decreasing rate of -6.5°C/km until the limit temperature -56.5°C. To validate that the captured energy would allow supplying a wireless sensor node, we performed the same experiment with a continuous current consumption of 220µA. It can be observed from Fig. 15 that sufficient energy is stored into the ultracapacitors to provide a regulated voltage over the whole flight, until t=8000s well beyond landing. This means that with appropriate power management, for example, periodic measurement and communication that would greatly reduce the power consumption, the same amount of energy should allow supplying the WSN node over a longer flight.

B. 0.8µm Smart Power SOI testchip

In this technology, we implemented the full battery-free multisource energy generator as shown in Fig. 16. Active diodes are also implemented in the voltage doubler used for the rectification of the signal originating from the piezoelectric transducer. The value of the external input capacitance of the voltage doubler, C in , is 10µF. This testchip, shown in Fig. 10, was designed with a regulated output voltage value Vout of 2.4V. The circuit was tested under the same conditions as the previous test chip. Firstly, we validated its functionality using one single energy source, the thermogeneration, using the thermal profile simulating a one-hour flight. The measurement data are presented in Fig. 17 that shows a voltage regulation at 2.45V operational at a time of 650s. To simulate the multisource energy generation, we replaced the piezoelectric transducer by its equivalent impedance (300Ω)

and applied an oscillating waveform with 60Hz frequency and 3V peak amplitude using a voltage generator. These parameters were defined in order to correspond to the experimental data of Fig. 3 in terms of frequency, and to a high level of vibrations (aircraft engines area and climb phase). As shown in Fig. 17, the regulated voltage is available much earlier (400s). However, this setup time is much higher than the one expected from the simulated one (a few seconds in Fig. 8). We analysed the reasons for this big discrepancy and found out that the ESD protection network, involving Zener diodes with a threshold of 6V, is unexpectedly triggered on and is deriving an important part of the energy generated by the piezoelectric transducer, thus considerably degrading the efficiency of the voltage doubler. To cope with this issue, the solution we contemplate consists in using a rectifier stage that is self-protected with regards to ESD without the addition of any protection devices. Despite this loss, the earlier start-up of the active diodes allows greatly improving the energy transfer to the ultracapacitors: up to 32% compared to the single TEG source. As a result the total energy collected is now 10J that corresponds to a power harvester efficiency, defined by (effective power on storage)/(harvester raw power on matched load), of 66%. The static efficiency of the standalone converter (without the ultracapacitors but with the small Cbias capacitor) was measured in two configurations (Fig. 18): by applying a DC power supply (3V or 4V) either to Vpol, (TEG and PIEZO source inputs being grounded) or directly to both TEG (Vin+/-) and PIEZO (Pz+/-) source inputs. The converter efficiency is computed by the ratio between the output and input powers measured using Agilent N6705B DC Power Analyzer. The remarkable result is that the active diode stages do not impact the global converter efficiency that is mainly defined by the LDO one. When the output voltage is stabilized, the quiescent current of this circuit is 300 nA for a DC power supply of 3.5V. Regarding the active diodes, they exhibited an average threshold of 26 mV. We also tested the circuit behaviour over temperature and obtained a regulated output voltage with a variation of 271 ppm/°C very similar to the one achieved with the 0.35µm CMOS technology.

Recently, several commercial solutions were proposed for a similar energy generation function. We compare in Table 1 the V. CONCLUSION

The multi-source energy generator presented in this paper exhibits two advantages compared to available commercial solutions. First, it allows using two sources of energy, which is a significant advantage for the concept of energy harvesting since the availability of a single energy source in a given environment is not guaranteed. The second advantage is its ultra-low power consumption that is critical for the autonomy of a battery-free wireless sensor node.

This circuit was designed for a specific SHM aeronautics application where transient gradients and permanent mechanical vibrations are used as energy sources. Vibrations provide a secondary source mainly used to self-bias the circuit, which does not include any battery and stores energy into ultracapacitors. The proposed circuit, realized in a 0.35µm CMOS technology, exhibits 200nA quiescent current and 300nA in the 0.8µm Smart power SOI technology well below values of commercial products. This system was also validated on a test bench emulated with realistic mission profiles. Finally, the innovative building blocks presented in this paper may find applications in field outside aeronautics.

Fig. 1 .

 1 Fig. 1. Cross-section (a) and photograph (b) of the experimental setup for the testing of the thermo electric generator: airflow is temperature-regulated and a small air volume in the capsule allows absorbing the water volume variation at phase change. View (b) shows the Thermostream® climatic chamber test head above the mockup.

Fig. 2 Fig. 2 :

 22 Fig.2shows a typical voltage response of a TEG under transient thermal gradients created with the apparatus of Fig.1including 10g of water and the climatic chamber described in section IV. This figure also provides the temperature profile applied to the TEG that corresponds to a typical one-hour flight. As can be observed, gradient polarity is reversed during descent, asking for rectifying the TEG output.

Fig. 3 :

 3 Fig. 3: Experimental amplitude spectrum of acceleration versus frequency. Acceleration was measured by an accelerometer located in aircraft engine area during descent. Both scales are logarithmic. The fundamental frequency component ranges typically between 40 and 60 Hz depending upon flight phase and engine model

Fig. 4 :Fig. 5 :

 45 Fig. 4: Electrical schematic of the proposed battery-free multisource energy generator connected to a sensor node. Diode symbols (Dp1, Dp2, Dn1, Dn2, Dp3) labelled with * are active diodes (see electrical schematic in Fig. 5) and capacitor symbols (UC1 & UC2) labelled with * are ultracapacitors. PTAT/BGR block provides the bias current for the active diodes and the voltage reference for the LDO. The energy transducers are represented by their respective electrical model.

Fig. 6 :

 6 Fig. 6: Symbols (top) and electrical schematic for the active diodes Dn (left) and Dp (right).

Fig. 7 :

 7 Fig. 7: Proposed nanowatt voltage reference.

Fig. 8 :

 8 Fig. 8: Comparative simulation of the generated voltage Vbias across the storage device Cbias during a one-hour flight. (a) Simulation set-up shown for the case of the combined energy sources, Pz input being grounded only the TEG source is used. The energy sources are simulated using their equivalent electrical model and a voltage input corresponding to the output of the respective transducer under typical flight conditions. (b) Simulation plot of Vbias voltage across

Fig. 9 :

 9 Fig. 9: Electrical schematic of the LDO with the proposed voltage reference circuit and the capability of generating two regulated output voltages, Vref or 2Vref, according to the command on transistor P6.

Fig. 10 :

 10 Fig. 10: Optical views of the two test chips: 0.35µm High-voltage CMOS (left) and 0.8µm smart power SOI (right).

Fig. 4 :

 4 Fig. 4: Electrical schematic of the ultra low-power (ULP) circuit with the external elements needed for the energy generation.

Fig. 12 :

 12 Fig. 12: Experimental characterization of the active diodes.

Fig. 13 :

 13 Fig. 13: Temperature behaviour of the regulated output voltage compared to the range defined by Monte Carlo simulations.

Fig. 14 presents

 14 Fig.14presents the results of measurements simulating a one-hour flight using the TEG for energy generation without drawing any current. As expected, since transient gradients build up at take-off, the circuit is only active after a few minutes (~10 min). The total energy that is collected during such a flight is equal to 4J. Compared to the original harvested energy of 15J obtained on a matched load, the losses are mainly due to the impedance mismatching and to the late triggering of the active diodes. The observed slight decrease of Vout voltage during cruise and on ground is related to the intrinsic leakage currents of

Fig. 14 :Fig. 15 :

 1415 Fig. 14: Energy generation with thermal gradient using a Micropelt MPG-D751 device: simulation of a one-hour flight. This test is performed without drawing current.

Fig. 16 :

 16 Fig. 16: Electrical schematic of the multisource and battery-free energy generator.

Fig. 17 :

 17 Fig. 17: Energy generation from thermal gradients using a Micropelt MPG-D751 device (_TEG) and from two sources, thermal gradients and vibrations (_TEG+PZ): experimental simulation of a one-hour flight. Vst_x is the voltage stored on the ultracapacitors and Vout_x is the output regulated voltage.

Fig. 18 :

 18 Fig. 18: Measured static converter efficiency using Agilent N6705B DC Power Analyzer for different output powers: comparison between setting the power supply, 3V or 4V, at the LDO source (ηLDO) and at the TEG/PIEZO inputs (η).

  performance of the proposed converter to the solution proposed by Linear Technology with the LTC3105 [25] and to the energy harvesting evaluation unit, TE-Power PLUS, proposed by Micropelt[26]. The results of research works[10] [13] [14] [27] are also given for reference. Although the different converters are difficult to compare in terms of complexity and performance, the quiescent current of the proposed converter is well below the value of both commercial and academic circuits. One of the circuits, namely the one from the Imperial College[START_REF] Kiziroglou | Heat storage power supply for wireless aircraft sensors[END_REF], has a very similar application is but only uses a thermogenerator for power harvesting. In terms of efficiency, at Vin=3V, despite the use of an LDO regulator, the converter we propose has a better quite honestly to the other circuits that use switched-mode DC/DC conversion.

Table 1 :

 1 Comparison of the performance of the proposed converter to commercial solutions LTC3105 from Linear Technology and TE-Power PLUS from Micropelt and to the results of research works [10] [13] [14] [27].

		LTC3105	TE-PowerPlus	[10]	[13] [27]	[14]	Proposed converter
	Quiescent current	10µA	3.2µA	NA	27µA	70µA	300nA
	Maximum
	efficiency Pout/Pin	90%	≥ 70%	NA	75-85%	50%	82%
	Harvested
	energy efficiency	NA	NA	60%	70%	NA	66%
	Regulated voltage	(LDO) 2.2V	2.4V	3.3V	1.2V	2.5V	2.4V

This work has been carried out within the framework of AUTOSENS and CREME projects funded by Fondation de Recherche