A new standardless approach in EPMA for quantitative analysis of uranium microparticles

Mouad Essani 1,2 Quantitative electron probe microanalysis (EPMA) is usually performed using analysis based on standards. The x-ray spectrum is acquired for an unknown sample and compared to the spectrum of a reference standard. Matrix effects are further corrected using conventional ZAF or 𝜑(𝜌𝑧) [START_REF] Philibert | L'analyse quantitative en microanalyse par sonde électronique[END_REF][2] methods. However, ZAF and 𝜑(𝜌𝑧) models were developed for bulk specimens with plane surfaces and involve many assumptions about the thickness and morphology of the sample. These models are no longer accurate when quantifying the composition of non-conventional samples such as microparticles. In fact, microparticles can have irregular shapes and exhibit high surface roughness and porosity. Thus, absorption path lengths and fluorescence inside the particle might be significantly different than those in a bulk polished sample. Moreover, quantitative models based on bulk standards do not take into account the loss of x-ray signal caused by electrons side scattering and/or transmission of electrons through the micro sized particles. Quantitative standardless procedures [START_REF] Trincavelli | Standardless quantification methods in electron probe microanalysis[END_REF] can be used to overcome these limitations. A virtual x-ray spectrum is calculated in that regard and used as a reference standard to quantify the sample composition. Standardless models based on Monte Carlo (MC) simulation were shown to be the most viable approaches to deal with non-conventional samples such as porous media [START_REF] Sorbier | EPMA of Porous Media: A Monte Carlo Approach[END_REF], thin films [START_REF] Hu | Method for the calculation of the chemical composition of a thin film by Monte Carlo simulation and electron probe microanalysis[END_REF] and microparticles [START_REF] Ritchie | Using DTSA-II to simulate and interpret Energy Dispersive spectra from Particles[END_REF].

In this work, we present a new MC quantitative standardless approach. Weight fractions are computed iteratively by means of the deviation between simulated and experimental intensities. The measured x-ray spectra are corrected after introducing the detector response function and compared to simulated ones (figure 3). The virtual x-ray spectra are calculated using the general-purpose MC code PENELOPE [START_REF] Salvat | A Code System for Monte Carlo Simulation of Electron and Photon Transport[END_REF] in which the morphological parameters of the sample can be described. PENELOPE uses Bote et al. ionization cross sections [START_REF] Bote | Calculations of inner-shell ionization by electron impact with the distorted-wave and plane-wave born approximations[END_REF] which was shown to give accurate results when M subshells are used [START_REF] Moy | Standardless quantification of heavy elements by electron probe microanalysis[END_REF].

A preliminary study was performed on a porous bulk sample of UO2 and a spherical microparticle of K411 (figure 1) glass analyzed at different beam incidence angles. The correction of intensity loss due to porosity, geometric (figure 2) and incidence angle effects is achieved, which leads to promising quantitative results with relative errors |𝐶𝑐-𝐶𝑛| 𝐶𝑛 <15%. The approach was further applied on micropaticles of U3O8, UF4 and UCeO.

Figure 1 Figure 2

 12 Figure 1 SEM image of an 8.4µm K411 glass spherical particle deposited on an aluminum substrate

Figure 3

 3 Figure 3 Comparison between measured and simulated spectrum for spherical particle of K411 glass. The deviation between intensities do not exceed 9%

  , Emmanuelle Brackx 1 , Philippe Jonnard 2 and Emmanuel Excoffier 1 1 CEA Marcoule, Laboratoire de Métallographie et d'Analyse Chimique, BP 17171, 30207 Bagnols sur Cèze, France

	2 Sorbonne Université, CNRS UMR 7614, Laboratoire de Chimie Physique-Matière et Rayonnement, 4 place Jussieu, F-75252 Paris
	Cedex 05, France


Gloria Faccanoni
email: gloria.faccanoni@univ-tln.fr

Cédric Galusinski
email: cedric.galusinski@univ-tln.fr

INFLUENCE OF A NONLINEAR DEGENERATE DIFFUSION ON AN ADVECTION-DIFFUSION EQUATION IN A DIFFUSE INTERFACE FRAMEWORK *

Keywords: phase change, Stefan problem, nonlinear degenerate diffusion

HAL is

Context

A heat exchanger is a device that facilitates the process of heat exchange between two fluids that are at different temperatures. Heat exchangers are used in many engineering applications, such as refrigeration, heating and air-conditioning systems, power plants, chemical processing systems, automobile radiators, and waste heat recovery units. Nuclear core of a pressurized water reactor is an example of heat exchangers.

In some applications, such as the flow in nuclear reactor cores, convection is characterized by a low Mach number, where the convective velocities are much slower than the speed of sound in the fluid. This has promoted the development of so-called low Mach number models, which filter the sound waves. In the context of pressurized water reactor cores, an asymptotic low Mach number model, called lmnc (for Low Mach Nuclear Core), has been derived and investigated in a series of papers [START_REF] Bernard | Study of a low Mach nuclear core model for two-phase flows with phase transition I: Stiffened gas law[END_REF][START_REF] Bernard | Study of a low Mach nuclear core model for single-phase flows[END_REF][START_REF] Dellacherie | On a low Mach nuclear core model[END_REF][START_REF] Dellacherie | Accurate steam-water equation of state for two-phase flow LMNC model with phase transition[END_REF][START_REF] Faccanoni | Thermal Diffusion and Phase Change in a Heat Exchanger[END_REF][START_REF] Faccanoni | Numerical Simulations of a low Mach number model in heat exchanger[END_REF][START_REF] Faccanoni | A homogeneous relaxation low Mach number model[END_REF]. The model was derived through an asymptotic expansion performed in the compressible Navier-Stokes equations with an energy source term. It consists of a transport equation upon a thermodynamic variable (here the total enthalpy), of a divergence constraint upon the velocity (with a nonlinear coupling source term which underlines the dilatation property of the flow) and of the momentum equation. The fluid is described by a single equation of state taking into account the phase transition by supposing that, when both vapour and liquid phases are present, they have the same pressure, temperature and chemical potential. The equation of state is piecewise defined with respect to the enthalpy identifying a liquid phase, a gas phase and a mixture of them. The thermal diffusion degenerates in the mixture region.

In the present paper we are interested in studying the influence of the thermal diffusion.

A heat exchanger model

For some bounded domain Ω ⊂ R 3 , the Navier-Stokes-Fourier system with source terms in a 3D conservative formulation reads

     ∂ t ρ + ∇ • (ρu) = 0, ∂ t (ρu) + ∇ • (ρu ⊗ u) + ∇p = ρg + ∇ • σ(u), ∂ t (ρh) + ∇ • (ρhu) = Φ + ∇ • (ω∇T ) + σ(u) : ∇u + ∂ t p + u • ∇p. (1)
The unknows are the velocity field (t, x) → u(t, x), the total enthalpy (t, x) → h(t, x) and the pressure (t, x) → p(t, x). The specific density ρ and the temperature T are related to the enthalpy h and the pressure p through an equation of state. The heat conductivity ω, assumed to be constant and isotropic for each phase, is also related to the enthalpy h and the pressure p. The power density (t, x) → Φ(t, x) ≥ 0 is a given function of time and space modelling the heating of the coolant fluid (e.g. due to the fission reactions in the nuclear core). Finally, g is the gravity field and σ(u) models viscous effects.

In a low Mach number regime (i.e. when the speed of fluid is much smaller than the speed of sound), an asymptotic expansions with respect to the Mach number can be applied to the compressible model, and this leads to a simpler model where the acoustic waves are filtered out. The "original" pressure p is decomposed in two pressure fields p(t, x) = p * + p(t, x), where p * is the reference state pressure (that is an average pressure constant in time and space) and p is the perturbational pressure (often called "dynamic pressure"). For low Mach number flows, an asymptotic analysis shows that p(t, x)/p(t, x) = O(M 2 ), where M is the Mach number. We can thus approximate p(t, x) by p * in the computation of thermodynamics quantities and we obtain the so-called LMNC model [START_REF] Bernard | Study of a low Mach nuclear core model for two-phase flows with phase transition I: Stiffened gas law[END_REF][START_REF] Dellacherie | On a low Mach nuclear core model[END_REF][START_REF] Faccanoni | Numerical Simulations of a low Mach number model in heat exchanger[END_REF], that can be written as

     ∂ t ρ + ∇ • (ρu) = 0, ∂ t (ρu) + ∇ • (ρu ⊗ u) + ∇p = ρg + ∇ • σ(u), ∂ t (ρh) + ∇ • (ρhu) = Φ + ∇ • (ω∇T ) + σ(u) : ∇u.
(2)

In this system the unknows are the velocity field u, the enthalpy h and the perturbational pressure (t, x) → p(t, x). The specific density ρ and the temperature T are linked, by an equation of state, to the enthalpy h and to the constant given thermodynamic pressure p * .

In the following the power density Φ is chosen constant and positive. We also neglect the viscous terms.

Diffusion term with phase transition

The fluid can be in liquid phase ℓ, vapour phase g or a mixture of them m.

if h ≤ h s ℓ (p), in vapor phase if h ≥ h s g (p)
and a mixture at saturation if h s ℓ (p) < h < h s g (p). Thus the global EoS is piecewise defined w.r.t. h s κ (p). In particular for the temperature we have (see Figure 2)

(h, p) → T (h, p) =      T ℓ (h, p), if h ≤ h s ℓ (p), T m (h, p) = T s (p), if h s ℓ (p) < h < h s g (p), T g (h, p), if h ≥ h s g (p).
Since the temperature is a function of the enthalpy h and the (thermodynamic) pressure p, the gradient of the temperature can be written as

∇T (h) = ∂T ∂h p ∇h + ∂T ∂p h ∇p.
In the lmnc model, the thermodynamic quantities are evaluated at p = p * , which is constant. Thus, from now letting drop the dependency upon p * , the diffusion term can be written as

ω(h)∇T (h) =      λ ℓ ∇h, if h ≤ h s ℓ , 0, if h s ℓ (p * ) < h < h s g , λ g ∇h, if h ≥ h s g ,
where λ κ def = ω κ /c p,κ for κ = ℓ or g with ω κ the heat conductivity (assumed to be constant and isotropic for each pure phase κ) and c p,κ def = ∂h ∂T p the isobar heat capacity of the pure phase κ. We note that if the open domain of the mixture is non-empty, a zero diffusion coefficient can be set in this region.

Let us denote

λ(h) def =      λ ℓ if h ≤ h s ℓ , λ m if h s ℓ < h < h s g , λ g if h ≥ h s g . ( 3
)
By defining λ m = 0 and by virtue of (2), with (3), the obtained equation on enthalpy comes down to

∂ t (ρh) -∇ • (λ(h)∇h) + ∇ • (ρhu) = Φ + σ(u) : ∇u.
It differs strongly on the enthalpy equation written in [START_REF] Voller | An implicit enthalpy solution for phase change problems: with application to a binary alloy solidification[END_REF] where the transition liquid-gas is forced at a sharp interface by the addition of a singular source term and a non-degeneracy of the thermal diffusion. The generality of the model proposed here is that the liquid-gas transition can also take place through a sharp interface or through a mixing region (in a diffuse interface framework). In the following, we will see that both can occur and when the liquid-gas transition is a sharp interface, our enthalpy solution coincides with that of the Stefan models written in temperature [START_REF] Rubinstein | The Stefan problem[END_REF] or enthalpy [START_REF] Voller | An implicit enthalpy solution for phase change problems: with application to a binary alloy solidification[END_REF]. This is the main point of this paper and an analytic solution, derived in the 1D steady-state case in the following section, shows this richness of the phase transition.

Steady-state model

In this section, we focus on the steady-state solution, i.e. we seek for h, u and p solution in

R 3 + of      ∇ • (ρu) = 0, ∇ • (ρu ⊗ u) + ∇p = ρg, ∇ • (ρhu) = Φ + ∇ • (λ(h)∇h).
In the unidimensional case, the first and last equations can be written as

∂ y (ρv) = 0, ∂ y (ρhv) = Φ + ∂ y (λ(h)∂ y h) in [0, +∞) (4)
with ρ a function of h. To model the configuration described in the Figure 1, the system is closed by the Dirichlet boundary condition corresponding to an inlet condition: the inlet enthalpy h(0) = h e < h s ℓ and the inlet flow rate (ρv)(0) = D e . Then the flow rate is constant on the domain:

(ρv)(y) = D e for all y ∈ R + .

(

In order to focus on a physically relevant solution, we assume that lim

y→∞ ∂ y (λ(h)∂ y h) = 0 so that, thanks to (4), lim y→∞ ∂ y h(y) = Φ D e .
We are thus interested in y → h(y) solution of the following non-linear strongly degenerate convectiondiffusion problem:

D e h ′ -(λ(h)h ′ ) ′ = Φ in R + (6) with
• Φ and D e two strictly positive constants,

• h → λ(h) the degenerate diffusion defined in (3),

• the two boundary conditions

h(0) = h e < h s ℓ [inlet Dirichlet boundary condition], lim y→∞ h ′ (y) = Φ De [asymptotic behaviour]. (7)
To better understand the role of each parameter, in the following we will consider three cases for the three constant values λ ℓ , λ m and λ g defined in (3): (3) λ m = 0 and λ ℓ , λ g > 0 (model with a degenerate diffusion) In the last case, which corresponds to the phase transition model, two cases have to be considered. Indeed, it will be seen that the mixture region can disappear for a sufficiently large gas diffusion coefficient λ g .

Steady solution without diffusion

Proposition 2.1 (Steady solution without diffusion). When λ ℓ = λ m = λ g = 0, the steady enthalpy solution of (6)-( 7) is

h(y) = h e + Φ D e y. ( 8
)
Proof. Without diffusion we have to solve D e h ′ (y) = Φ with h(0) = h e < h s ℓ . The solution is thus h(y) = h e + yΦ/D e . □

Let us denote x s ℓ (resp. x s g ) the liquid/mixture transition (resp. mixture/gas) without diffusion (see Figure 3). The fluid is in pure liquid phase for y ≤ x s ℓ , in pure vapour phase for y ≥ x s g and a mixture at saturation if x s ℓ < y < x s g where

x s ℓ = h s ℓ -h e Φ De , x s g = h s g -h e Φ De . (9)

Steady solution with a strictly positive diffusion

Without degeneracy of diffusion, the unique solution can be explicited as follows.

Proposition 2.2 (Steady solution with a positive diffusion). If λ(h), defined in (3), is piecewise positive constant for all h, the steady enthalpy solution of (6)-( 7) is

h(y) =                  h ℓ (y) def = C 1 + Φ D e y + C 2 exp y λ ℓ /D e if y ≤ y s ℓ h m (y) def = C 3 + Φ D e y + C 4 exp y λ m /D e if y s ℓ ≤ y < y s g h g (y) def = h s g + Φ D e (y -y s g ) if y ≥ y s g ( 10
)
where the constants C 1 , C 2 , C 3 , C 4 , depending on y s ℓ and y s g , are

C 1 = h e L s ℓ + Φ De (y s ℓ -1) L s ℓ -1 , C 2 = h e -C 1 , L s ℓ = exp y s ℓ λ ℓ /D e C 3 = Φ De y s ℓ -h s ℓ M s g -Φ De y s g -h s g M s ℓ M s ℓ -M s g , C 4 = Φ De y s g -h s g -Φ De y s ℓ -h s ℓ M s ℓ -M s g , M s κ = exp y s κ λ m /D e ,
for κ = ℓ or g and y s ℓ and y s g satisfies the continuity flux

λ ℓ h ′ ℓ (y s ℓ ) = λ m h ′ m (y s ℓ ), λ m h ′ m (y s g ) = λ g h ′ g (y s g ). ( 11
)
Proof. The non-degenerate diffusion, i.e. λ(h) > 0, leads to smooth solutions in the sense that discontinuities can not occur since order one distributions generated by the diffusive term and discontinuous solutions can not be compensated by order zero distribution. Furthermore, for continuous solutions, the transport term is then a diffuse measure leading to a continuous diffusion flux. Then, the continuous solution induces a mixture region and satisfies the flux continuity [START_REF] Myers | The Stefan problem with variable thermophysical properties and phase change temperature[END_REF].

Due to the growth of the steady state solution (Φ and D e are strictly positive), three regions, liquid, mixture, gas, are juxtaposed from low to high y values. On each region, an explicit solution solves the equation ( 6), [START_REF] Rubinstein | The Stefan problem[END_REF] where the unknows C 1 , C 2 , C 3 , C 4 , C 5 , C 6 , y s ℓ , y s g have to satisfy • the enthalpy continuity conditions at inlet and transition points:

h(y) =              h ℓ (y) def = C 1 + Φ De y + C 2 exp y λ ℓ /De if y ≤ y s ℓ h m (y) def = C 3 + Φ De y + C 4 exp y λm/De if y s ℓ ≤ y < y s g h g (y) def = C 5 + Φ De y + C 6 exp y λg/De if y ≥ y s g
h ℓ (0) = h e , h ℓ (y s ℓ ) = h s ℓ , ( 13
)
h m (y s ℓ ) = h s ℓ , h m (y s g ) = h s g , ( 14
)
h g (y s g ) = h s g , (15)
• the flux continuity conditions at transition points:

λ ℓ h ′ ℓ (y s ℓ ) = λ m h ′ m (y s ℓ ), λ m h ′ m (y s g ) = λ g h ′ g (y s g ), (16)
• the asymptotic behaviour

lim y→∞ h ′ g (y) = Φ D e . ( 17
)
The conditions ( 13), ( 14), ( 15) and (17) allow to define C 1 to C 6 as functions of the transition points y s ℓ and y s g . These points are thus implicitly defined by ( 16) and finally we obtain the solution. We note that the transparent boundary condition (17) implies C 6 = 0 so that y → h g (y) is linear and compatible with the semi-infinite domain. Moreover, although y → h ′ g (y) does not involve λ g , the flux continuity (16) does. □ The Figure 4 shows the graph of y → h(y) with λ m > 0, that is a non-degenerate artificial diffusion in the mixture region, leading to a continuous solution. Depending on the increasing ratio λ g /D e when other parameters are fixed, the mixture region reduces and is pushed on the left. A stiffer solution appears at the gas transition. On the last plot (Figure 4c), the mixing region is very sharp and is not included in the large mixing region occurring without diffusion in the whole phases.

The first plot of the Figure 5 is the one of the Figure 4a. In the next figures is plotted a reduced diffusion coefficient λ m . It is then exhibited the convergence of the enthalpy as λ m goes to zero. A discontinuous solution is suggested at the limit. It occurs at the transition to the gas phase. For such diffusions in the liquid and gas phases, the mixture region is slightly affected by the decreasing diffusion λ m and remains at the limit.

We will see in the following section that, when the diffusion is allowed to degenerate, i.e. the diffusive coefficient λ(h) vanishes for values of h in the mixture region, solutions are not necessarily smooth and discontinuous weak solutions must be sought.

Steady solution with degenerate diffusion

The degenerate diffusion can lead to a jump in the diffusion flux at phase transition. For example, as the enthalpy flux is zero in the mixture region, a flux diffusion jump occurs at the interface with the mixture if the space derivative of the enthalpy is not zero at this transition in the liquid (or gas) phase. A Dirac measure follows and has to be compensated by an other Dirac measure in the transport term. A discontinuous solution is then necessary. Due to the growth of the steady state solution (Φ and D e are positive) and to the presence of the liquid phase at inlet, a liquid phase is present for low y values, then follows a transition which can be liquid-mixture (followed by the transition mixture-gas) or liquid-gas. The two possible cases are presented below.

Case I: Liquid/mixture/gas

The solution with mixture is possible for a sufficiently low ratio λg De Φ De and is explicited in the following proposition.

Proposition 2.3 (Steady solution with degenerate diffusion: liquid/mixture/gas case). For all λ ℓ > 0, if

λ g D e Φ D e < h s g -h s ℓ ( 18
)
the mixture zone is present and the solution of (6)-( 7) satisfies

h(y) =                  h ℓ (y) def = h e + Φ D e y + Φ D e λ ℓ D e 1 -exp y λ ℓ /D e exp -y s ℓ λ ℓ /D e if y ≤ y s ℓ h m (y) def = h s ℓ + Φ D e (y -y s ℓ ) if y s ℓ ≤ y < y s g h g (y) def = h s g + Φ D e (y -y s g ) if y ≥ y s g ( 19
)
where

• the position y s ℓ separating the liquid phase to the mixture phase is implicitly defined by h ℓ (y s ℓ ) = h s ℓ ; • the position y s g separating the mixture phase and the gas phase is computed w.r. Proof. As for the nondegenrate case, on each phase or mixture, the integration of the solution gives

h(y) =                  h ℓ (y) def = C 1 + Φ D e y + C 2 exp y λ ℓ /D e if y ≤ y s ℓ , h m (y) def = C 3 + Φ D e (y -y s ℓ ) if y s ℓ ≤ y < y s g , h g (y) def = C 4 + Φ D e (y -y s g ) + C 5 exp y λ g /D e if y ≥ y s g . ( 20
)
The five constants C 1 • • • C 5 have to be determine and a particular choice will be made to define a unique viscosity solution. Searching a solution of ( 6) with mixture, the transition liquid-mixture induces a free dissipative flux in the mixture. A positive jump on dissipative flux can not be compensated by the enthalpy jump on the transport term since the enthalpy is increasing. As a matter of fact, if we denote h the enthalpy jump at the liquid-mixture interface, the jump relation of ( 6) leads to

D e h -0 × h ′ (y s,+ ℓ ) + λ ℓ h ′ (y s,- ℓ ) = 0, that is λ ℓ h ′ (y s,- ℓ ) = -D e h .
Since y → h(y) has to be increasing, we should have h ≤ 0. Hence, in the present case, the only possibility is to have both h = 0 and h ′ (y s,- ℓ ) = 0. The inlet Dirichlet boundary condition and the Dirichlet condition h(y s ℓ ) = h s ℓ define the constants C 1 and C 2 parametrized with respect to y s ℓ . Then, the additional Neumann boundary condition at y s ℓ defines implicitly the unique value of y s ℓ , C 1 and C 2 , that is the unique solution in the known liquid phase.

The expressions of the solution in the mixture region follow by the continuity of the enthalpy at the interface with liquid, thus C 3 = h s ℓ . Furthermore, since C 5 = 0 by virtue of the asymptotic behavior [START_REF] Faccanoni | Thermal Diffusion and Phase Change in a Heat Exchanger[END_REF],

h g (y) = C 4 + Φ D e (y -y s g ).
The difficulty is then to define the transition point y s g and the constant C 4 . The jump relation at this point is obtained in the same way as for the point y s ℓ , but this time, the two jumps can compensate each other:

D e h -λ g h ′ (y s,+ g ) + 0 × h ′ (y s,- g ) = 0, that is h = λ g D e Φ D e ,
which has to be lower than the maximal gap h s g -h s ℓ . There is an infinity of solutions to such a problem in the sense of distribution. The transition y s g separating enthalpy lower than h s g to upper can cover an interval, that is,

h g (y s g ) ∈ h s g , h s g + λ g D e Φ D e .
Thus,

h m (y s g ) = h g (y s g ) -h = h g (y s g ) - λ g D e Φ D e .
In order to define an unique solution, we choose a solution corresponding to the smallest mixture region, that is with the smallest value of enthalpy h g (y s g ) = C 4 = h s g in the gas region at the transition. The reader can convince himself (see Figure 5) that this choice corresponds to the viscosity solution by studying the limit of the viscous solution obtained in Proposition 2.2. This implies that

h m (y s g ) = h s g - λ g D e Φ D e
and the position y s g follows

y s g = y s ℓ + D e Φ (h s g -h s ℓ ) - λ g D e .
□ Remark 2.4. Gas diffusion reduces the mixture region for steady solution since

(y s g -y s ℓ ) = (x s g -x s ℓ ) - λ g D e
where x s κ are the transition points without diffusion defined in [START_REF] Faccanoni | A homogeneous relaxation low Mach number model[END_REF]. We remark also that y s ℓ > x s ℓ if λ ℓ > 0. We point out that if (and only if) the condition

λ g D e Φ D e ≤ h s g -h s ℓ
is not fulfilled, a transition liquid-gas has to be written (cf. the following section).

The Figure 6 shows some steady solutions when (18) is fulfilled so that the mixture zone is present. The mixture region is always smaller than the one obtained without diffusion, i.e. y s g -y s ℓ ≤ x s g -x s ℓ . More precisely, we observe:

• The Figure 6a shows the steady solution when Φ = D e , λ ℓ = 2D e and λ g = 0. We can see the effect of the diffusion only in the liquid phase: the liquid/mixture transition occurs at y s ℓ > x s ℓ and the slope of h ℓ at y s ℓ is zero. • Figure 6b shows the steady solution when Φ = D e , λ ℓ = 0 and λ g = D e . We can see the effect of the diffusion only in the vapour phase: the mixture/vapour transition occurs at y s g < x s g ; the width of the mixture zone is reduced w.r.t. the case without diffusion and there is a jump on the enthalpy h at transition y s g . • Finally Figure 6c shows the steady solution when Φ = D e , 2λ ℓ = D e and λ g = D e . We can see the effect of the diffusion both in the liquid and in the vapour phases: the liquid/mixture transition occurs at y s ℓ > x s ℓ and the slope of h ℓ at y s ℓ is zero; the mixture/vapour transition occurs at y s g < x s g ; the width of the mixture zone is reduced w.r.t. the case without diffusion as observed in the remark 2.4 and there is a jump on the enthalpy h at transition y s g .

Case II: Liquid/gas (no mixture, Stefan problem)

The solution without mixture occurs for a sufficiently large ratio λg De

Φ

De and is explicited in the following proposition.

Proposition 2.5 (Steady solution with degenerate diffusion: liquid/gas case). For all λ ℓ > 0, if the mixture zone is not present and the solution of (6)- [START_REF] Faccanoni | Thermal Diffusion and Phase Change in a Heat Exchanger[END_REF] satisfies

λ g D e Φ D e ≥ h s g -h s ℓ ( 21
)
h(y) =        h ℓ (y) def = h e + Φ D e y + (h s g -h s ℓ ) - λ g D e - λ ℓ D e Φ D e 1 -exp y λ ℓ /D e exp -y s λ ℓ /D e if y < y s h g (y) = h s g + Φ D e (y -y s ) if y > y s
where the position y s separating the liquid phase to the gas phase is implicitly defined by h ℓ (y s ) = h s ℓ . Proof. As shown in the proof of Proposition 2.3, there isn't viscosity solutions to ( 6)-( 7) with mixture if

λ g D e Φ D e ≥ h s g -h s ℓ .
On each phase, the integration of the differential equation gives

h(y) =          h ℓ (y) def = C 1 + Φ D e y + C 2 1 -exp y λ ℓ /D e if y < y s h g (y) def = C 3 + Φ D e (y -y s ) + C 4 exp y λ g /D e if y > y s (22)
The four constants C 1 . . . C 4 and the transition point y s have to be determined. The inlet Dirichlet boundary condition and the asymptotic condition (7) define the constants C 1 = h e and C 4 = 0.

A particular choice will be made to define an unique viscosity solution. Such a viscosity solution has to satisfy h(y s,-) = h s ℓ , h(y s,+ ) = h s g . This implies that C 3 = h s g and

C 2 = h s ℓ -h e -Φ De y s 1 -exp y s λ ℓ /De .
The jump relation at the transition point y s reads

D e (h s g -h s ℓ ) = λ g h ′ (y s,+ ) -λ ℓ h ′ (y s,-). Since C 4 = 0, h ′ g (y) = Φ De for all y so that λ ℓ D e h ′ ℓ (y s,-) = λ g D e Φ D e -(h s g -h s ℓ ),
which is positive under the assumption (21) and is then compatible with the increasing enthalpy h.

Computing the derivative of the enthalpy in the liquid phase, we obtain an implicit definition of the position y s :

(h s g -h s ℓ ) + Φ D e λ ℓ D e - λ g D e exp y s λ ℓ /D e -1 = (h e -h s ℓ ) + Φ D e y s exp y s λ ℓ /D e .
We now prove that this value is unique. Let us we denote t def = exp

y s λ ℓ /De ≥ 1, C def =(h s g -h s ℓ ) + Φ De λ ℓ De - λg De and R = Φ De λ ℓ
De . The equation on the flux continuity can be rewritten as C(t -1) = [(h e -h s ℓ ) + R ln(t)] t, thus the position y s is implicitly defined by the zero of the function

f (t) = (C + h s ℓ -h e )t -Rt ln(t) -C, t ≥ 1.
Since f (1) = h s ℓ -h e > 0, f (t) ---→ t→∞ -∞, and f ′′ (t) = -R/t < 0, the function f admits an unique zero for t > 1. □

Figure 7 shows some steady solutions when λg De Φ De ≥ h s g -h s ℓ so that the mixture zone is not present. • Figure 7a shows the steady solution when Φ = D e , λ ℓ = D e and λ g = 2D e . This corresponds to the smallest choice of λ g such that the jump is equals to h s g -h s ℓ . The mixture is not present, the liquid/vapour transition occurs at y s and the slope of h ℓ at y s ℓ is zero. • Figure 7b shows the steady solution when Φ = D e , λ ℓ = D e and λ g = 4D e . The jump is equals to h s g -h s ℓ , the mixture is not present, the liquid/vapour transition occurs at y s and the slope of h ℓ at y s ℓ is strictly positive. When λ ℓ → 0 the slope tends to +∞:

h ′ ℓ (y s,-) = λg De Φ De -(h s g -h s ℓ ) λ ℓ De ≥ 0 and ---→ λ ℓ →0

+∞

• Finally Figure 7c shows the steady solution when 2Φ = D e , λ ℓ = 0 and λ g = 5D e . The mixture is not present, the liquid/vapour transition occurs at y s . Since λ ℓ = 0 there is a jump greater to h s g -h s ℓ so that we have a jump in liquid phase that corresponds to a jump on temperature and h ℓ (y s ) < h s ℓ :

h (y s ) = h s g -h ℓ (y s ) = λ g D e Φ D e ≥ h s g -h s ℓ .

Generalization to negative power densities

An enthalpy jump was exhibited at the gas transition point y s g . The gas phase does not play a different role from the liquid one: it is the situation of a heating a fluid (that is Φ > 0) that induces an increasing solution in the direction of flow with a jump localized at the last transition point.

In the case of a cooling heat exchanger (that is Φ < 0), the fluid is injected in gas phase (h e > h s g ): the stationary solution is thus decreasing in the direction of flow and it presents:

• a jump at the mixture-liquid transition point in case of

λ g D e |Φ| D e < h s g -h s ℓ ;
• a discontinuity at the gas-liquid transition point otherwise.

Link to the Stefan problems

The solution exhibited in the Proposition 2.5 is exactly the stationary solution of the following stationary Stefan problem, a free boundary problem. That is, find the position y s so that the over-determinate elliptic problem is satisfied: Written with the temperature as unknown, we obtain the classic stationary Stefan problem, with a continuous temperature at the transition and a jump condition on thermal flux, as it can be found in [START_REF] Rubinstein | The Stefan problem[END_REF],

                     D e h ′ -λ ℓ h ′′ = Φ in ]0, y s [ h(0) = h e < h s ℓ , h(y s,-) = h s ℓ , D e h ′ -λ g h ′′ = Φ in ]y s , +∞[ h(y s,+ ) = h s ℓ , lim y→∞ h ′ (y) = Φ De , λ g h ′ (y s,+ ) -λ ℓ h ′ (y s,-) = D e (h s g -h s ℓ ). (23)
                     c p,ℓ D e T ′ -ω ℓ T ′′ = Φ in ]0, y s [ T (0) = T e < T s , T (y s,-) = T s , c p,g D e T ′ -ω g T ′′ = Φ in ]y s , +∞[ T (y s,+ ) = T s , lim y→∞ T ′ (y) = Φ cp,gDe , ω g T ′ (y s,+ ) -ω ℓ T ′ (y s,-) = D e (h s g -h s ℓ ). (24)
Note that the last relation of (24) describes the jump on heat flux due to latent heat. This problem can be solved even if (18) holds, but the solution would loose a physical meaning due to inconsistent temperature with respect to the considered phase. As a matter of fact, the last relation of ( 24) gives

ω ℓ T ′ (y s,-) = λ g Φ D e
-D e (h s g -h s ℓ ) < 0 under (18), which induces a temperature above T s in a neighborhood of y s in the liquid phase.

To summarize, the Stefan formulations (23) or (24) have to be used only under (21) contrary to ( 6)-( 7) which is valid under (21) and (18). Under (21), the solutions of the formulations are the same; under (18), the Stefan solution of (23) or (24) isn't physical. That is why the PDE formulation ( 6)-( 7) has to be preferred for numerical solving as it is shown in Appendix A, extended to the non-stationary case.

Conclusion and perspectives

It has been shown the expression of the physical stationary solution for a one dimensional degenerate nonlinear elliptic problem modelling enthalpy subjected to a flow and thermal diffusion with phase change. The expression of the solution depends strongly on the ratio of heat deposit with flow rate to compare with the ratio of flow rate and diffusion coefficient in the gas phase (see (21)). It is notable that discontinuous solutions appear in the context of a modelization with diffuse interface. Furthermore, this mixture region which always exists without thermal diffusion, can disappear depending on the ratio (21). The selected solution is the viscosity solution leading to a minimal mixture region. The higher the gas diffusion coefficient, the smaller the mixing region, until it disappears. Some unsteady simulations are proposed in the appendix but a more thorough numerical study is to follow for the full model proposed (2) in a future work.

A. Numerical Time-Dependent solution of the constant-density model

The time-dependent 1D model ( 4) can be written in a non-conservative formulation as

∂ t ρ + ∂ y (ρv) = 0, ρ∂ t h + ρv∂ y (h) = Φ + ∂ y (λ(h)∂ y h), ∈ R + × R + .
In the stationary case, the model does not depend on the expression of the density with respect to the enthalpy. It is remarkable that whatever the density expression, the time-asymptotic solution depends only on the uniform flow rate. In the physical model, the density is assumed to be a continuous function of the enthalpy h, piecewise defined for each phase or mixture. In the transient regime, the law describing the density with respect to the enthalpy affects the unsteady solution. With the EoS describing a mixture at saturation as well as the pure phases, the density is a continuous and piecewise defined function of the enthalpy (and the constant

• Finally Figure 9c shows the particular case of a zero diffusion in the liquid phase leading to an increasing jump at the transition point to the gas phase. For sufficiently large time, the jump exceeds the value h s g -h s ℓ due to the lack of diffusion in the liquid phase which behaves like the mixture.

Figure 1 .

 1 Figure 1. The heat exchanger configuration

Figure 2 .

 2 Figure 2. Temperature T as a function of the enthalpy h at a constant pressure p *

Figure 3 .

 3 Figure 3. Steady solution without diffusion

  (a) λg = De (b) λg = 2De (c) λg = 4De

Figure 4 .

 4 Figure 4. Non degenerate diffusion for all h. Influence of λ g . Common parameters: Φ = D e , 2λ ℓ = D e , 10λ m = D e

Figure 5 .

 5 Figure 5. Non degenerate diffusion for all h. Convergence of λ m → 0. Common parameters: Φ = D e , 2λ ℓ = D e , λ g = D e

( a )Figure 6 .

 a6 Figure 6. Degenerate diffusion. Case I: Mixture zone at steady state. Common parameters: Φ = D e and λ m = 0.

Figure 7 .

 7 Figure 7. Degenerate diffusion. Case II: No mixture zone at steady state. Common parameter: λ m = 0.

  solution (a) Diffusion only in the liquid phase: λ ℓ = 2De, λg = 0. We observe the influence on y s ℓ and the slope at y s ℓ . solution (b) Diffusion only in the gas phase: λ ℓ = 0, λg = De. We observe the influence on ∆y s κ the width of the mixture zone and jump at y s g . solution (c) Diffusion both in liquid and vapour phases: λ ℓ = De, λg = De.We observe the existence of the mixture zone at steady state.

Figure 8 .

 8 Figure 8. Mixture zone at steady state. Common parameters: Φ = D e , λ m = 0.

  solution (a) Φ = De, λ ℓ = De, λg = 2De. We observe the disappearing of the mixture zone at steady state solution (b) Φ = De, λ ℓ = De, λg = 4De. We observe the disappearing of the mixture zone and the slope at y s,- can be > Φ solution (c) 2Φ = De, λ ℓ = 0, λg = 5De. Diffusion only in vapour phase, no mixture zone, jump in liquid phase at y s

Figure 9 .

 9 Figure 9. Diffusion both in liquid and vapour phases. No mixture zone at steady state.

This work has been achieved within the framework of the M2SIR project, financially supported by CNRS AAP "NEEDS 2020", and its partners Andra, BRGM,

thermodynamic pressure p * ). For simplicity, we choose here to exhibit the transient regime in the case of an artificially uniform in space and constant in time density. Let us consider ρ(h) = r > 0 constant for all h. The system becomes

with the boundary condition (rv)(0) = D e > 0 and h(0) = h e < h s ℓ . Thus v(y) = v(0) for all y ∈ R + and (ρv)(y) = rv(0) = D e for all y ∈ R + . The enthalpy equation can be written as

The unsteady model, written on the domain (t, x) ∈]0, T [×]0, +∞[, reads

, associated to the following boundary and initial conditions,

We define δt the time step and h n ≈ h(nδt). It is then essential to numerically solve this evolution problem by a fully implicit scheme,

associated to a gradient scheme [START_REF] Eymard | Gradient schemes for the Stefan problem[END_REF], in order to produce a viscosity solution. All the simulations respects the enthalpy properties at transition (jump relations), namely a jump of enthalpy (when it has to be) at the transition to the gas phase with an enthalpy close to h s g (up to the order of space discretization step) at the first node where enthalpy is upper than h s g . The increasing solution is also observed on the discrete solution avoiding spurious numerical oscillations.

Note that the solution converges, when time goes to infinity, to the analytic solution explicited in the Propositions 2.2 and 2.3. Note also that mixture appears in transient regime even if the asymptotic solution solves the Stefan problem with a transition liquid/gas, as on the Figure 9.

On Figures 8 and9, the initial condition is in liquid phase, then appears the mixture and finally the gas phase. The parameters of Figure 8 satisfy condition (18) with mixture on time-asymptotic solution.

• Figure 8a shows the solution without diffusion in the gas phase and the solution shows the characteristic transition liquid-mixture with a zero space derivative in the liquid at the motionless transition. • Figure 8b shows the solution without diffusion in the liquid phase and the solution shows the characteristic transition mixture-gas with a growing jump of enthalpy at the moving transition. • Finally Figure 8c shows the solution with diffusion in the liquid and gas phase gathering the two previous behaviors. The parameters of Figure 9 satisfy condition (21) without mixture on the time-asymptotic solution.

• Figure 9a shows the apparition of mixture then gas phase. The jump on the enthalpy is increasing along time and is asymptotically the maximal jump h s g -h s ℓ . This is the critical case where mixture disappears asymptotically with a zero space derivative in the liquid at the transition.

• Figure 9b shows the same phenomena with a bigger diffusion in the gas phase leading to positive slope in the liquid at the transition for sufficient large times.


Philippe Jonnard
email: philippe.jonnard@sorbonne-universite.fr

Jiaping Zhang

Christopher Pendenque

Karine Le Guen

Renaud Delaunay

Ian Vickridge

Didier Schmaus

Jiaping Zhang

Etude De Guides D'ondes

J.-P Zhang

Etude de guides d'ondes par PIXE en mode Kossel

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

L'émission de rayons X induite par protons (PIXE) combinée à la diffraction de Kossel est développée comme méthode de caractérisation structurale non-destructive de couches minces d'épaisseur nanométrique. Cette méthode consiste à mesurer en fonction de l'angle de détection (grazing exit) l'intensité d'une ou plusieurs émissions X caractéristiques d'éléments présents dans un empilement de couches minces. Ici, l'empilement est choisi de sorte qu'il se comporte comme un guide d'onde planaire, c'est-à-dire que dans une certaine gamme angulaire des interférences constructives se développent, conduisant à la variation en profondeur l'intensité du champ électrique associée à la radiation émise. La mesure PIXE-Kossel menée autour de la gamme angulaire du guide d'onde permet d'obtenir une information sensible en profondeur et ainsi de remonter à la structure l'empilement en particulier de ses interfaces.

La méthode PIXE-Kossel est reliée aux méthodes basées sur l'observation d'intensités d'émissions en mode d'onde stationnaire dans des empilements périodiques [START_REF] Le Guen | Kossel effect in periodic multilayers[END_REF]. Nous montrons ici qu'elle peut se généraliser à des systèmes non-périodiques de couches minces, pourvu qu'un système de guide d'ondes puisse être créé à l'intérieur de la structure. Nous montrons aussi que l'emploi d'une caméra à rayons X à dispersion d'énergie permet de mesurer l'intensité des émissions X en fonction de l'angle de détection, que nous appelons courbe de Kossel, avec un dispositif expérimental simple ne nécessitant aucun mouvement mécanique, Fig. 1.

La méthode est appliquée à l'étude d'un empilement Pt/Fe/Pt déposé par pulvérisation magnétron sur un substrat de silicium [START_REF]Kossel interferences of proton-induced X-ray emission lines to study thin film waveguides[END_REF]. Les épaisseurs visées des couches étaient comprises entre 4 et 10 nm et choisies de sorte qu'un système de guide d'onde puisse se développer aux énergies de photon des émissions caractéristiques Fe K et Pt L, Fig. 2. Un faisceau de protons de 1,6 MeV généré par l'accélérateur Van de Graaff de la plateforme SAFIR (Système d'analyse par Faisceaux d'ions Rapides) a été utilisé pour produire l'ionisation des atomes de l'échantillon. Une caméra CCD ANDOR iKon-M à dispersion d'énergie, 1024 × 1024 pixels de 13 × 13 µm, a été utilisée pour enregistrer les émissions de rayons X générées. Cette caméra CCD produit un spectre de rayons X en chacun de ses pixels. Dans notre cas, chaque colonne de pixels correspond à un angle de détection, Fig. 1.

Nous montrons sur la figure 3 les courbes de Kossel obtenues aux énergies des émissions Fe K et Pt L ainsi que leurs simulations. De fortes variations des intensités sont observées dans la gamme dans la gamme angulaire où le guidage d'onde prend place, avec un maximum vers 0.6° et 0.4° aux énergies des émissions Fe K et Pt L respectivement, en accord avec les distributions d'intensité du champ électrique montrées en Fig. 2. Si l'empilement conçu initialement, une tricouche, est choisi pour simuler les courbes, l'accord est moins bon que si un composé Pt69Fe31 épais d'environ 1,5 nm est considéré aux interfaces Pt-sur-Fe et Fe-sur-Pt. Un tel composé est attendu dans le système Pt-Fe [START_REF] Sato | Magneto-optical spectra in Pt/Co and Pt/Fe multilayers[END_REF]. Les mesures de réflectivité X, sensibles à la structure de l'empilement, ne donnent pas de différence claire entre les types de structures, tri-couche avec ou sans composé interfacial. L'ajout de l'information élémentaire, apportée par PIXE-Kossel, permet d'apporter une description réaliste de l'empilement de couches minces d'épaisseur nanométrique.

Figure 1 .

 1 Figure 1. Schema d'une expérience PIXE-Kossel.

Figure 2 .

 2 Figure 2. Simulation de la distribution angulaire et en profondeur du champ électrique dans un empilement Pt/Fe/Pt aux énergies des émissions (a) Fe Kα (6403 eV) et (b) Pt Lα (9442 eV). La surintensité dans la couche de fer démontre l'effet du guidage d'onde à 0.6° (a) ou 0.4° (b).

Figure 3 .

 3 Figure 3. Comparaison des corubes de Kossel expérimentales (lignes bleues) et simulées aux énergies des émissions Fe Kα (gauche) et Pt Lα (droite) dans l'empilement Pt/Fe/Pt montré en Fig. 1. Les simulations sont effectuées soit avec les paramètres (épaisseur et nombre de couches) de l'empilement définis lors de la conception du guide d'onde (courbes noires), soit en considérant un composé interfacial FePt (courbes rouges).


Zeaiter Hussein
email: hussein.zeaiter@femto-st.fr

Baala Oumaya
email: oumaya.baala@femto-st.fr

Francois Spies
email: francois.spies@femto-st.fr

Val Thierry
email: thierry.val@univ-tlse2.fr

Measuring accurate Angle of Arrival of weak LoRa signals for Indoor Positionning

Keywords: Autocorrelation function (ACF), peak detection, Fast Fourier Transform (FFT), Software Defined Radio (SDR), Angle of Arrival (AoA), Long-Range (LoRa), Negative SNR USRP Source Sample rate: 2M Centre frequency: 868M Bandwidth: 1M 2.048k FFT FFT size: 2.048k Window: Blackmanharris Forward/Reverse: Forward Multiply Conjugate Vec Length: 2.048k FFT FFT size: 2.048k Window: Blackmanharris Forward/Reverse: Reverse Complex to Mag Vec Length: 2.048k 2.048k Max Threshold Low: α High: α Initial State: 0 FFT FFT size: 2.048k Window: Blackmanharris Forward/Reverse: Forward

In this paper, we propose an Autocorrelation method for measuring the angle of arrival (AoA) of a weak LoRa signal. A weak LoRa signal has a negative SNR down to -20 dB. The objective is to detect a LoRa signal that operates at low transmission power (TX). Operating at low transmission power (TX) reduces power consumption and extends the battery life of LoRa devices. Besides, the transmission of weak signals strengthens the radio communication protocol, preventing an enemy device from accessing the location coordinates. The detecting algorithm consists of finding Autocorrelation peaks of the LoRa signal. We show that Autocorrelation peaks decrease when the signal is buried in the noise. However, using a large number of Fast Fourier Transform (FFT) will increase the Autocorrelation peaks and the signal-to-noise ratio (SNR). Once the peak of the LoRa signal is detected under the noise, the algorithm will calculate the AoA. All of the proposed algorithms are implemented using a Universal Software Radio Peripheral (USRP), Software Defined Radio (SDR) receiver with the help of GNU Radio software. We, therefore, believe that our Autocorrelation method can detect the LoRa signal accurately and measure the AoA at very low SNR in real-time, being usable for indoor positionning.

I. INTRODUCTION

Nowadays, The Internet of Things (IoT) is experiencing rapid growth. The IoT is a network of connected devices that enables connectivity and data exchange for a wide range of applications, including smart cities, healthcare, industrial and infrastructural applications, transportation, and navigation [START_REF] Hassan | Internet of Things and Its Applications: A Comprehensive Survey[END_REF]. IoT applications face many challenges, including security, low-cost devices, long-range communication, low data rate, This work is supported by the Femto-st institute in collaboration with Jolliet-Innovation and SECA companies in the framework of the SYLOIN FEDER BPI Project. and low power consumption. Today, we can address these challenges with the new low power wide area network (LP-WAN) paradigm. LPWAN standards enable wide communication coverage, low data rates, and long battery life. These advantages led the researchers to investigate the performance of the LPWAN technology both outdoors and indoors [START_REF] Muzaffar | Longrange wireless sensor networks for geo-location tracking: Design and evaluation[END_REF][START_REF] Zafari | A survey of indoor localization systems and technologies[END_REF]. Many LPWAN technologies have been studied recently, such as LoRaWAN, Sigfox, and NB-IoT. LoRaWAN is one of the latest LPWAN technologies that has attracted many researchers interest in outdoor and indoor navigation [START_REF] Choi | Low-Power LoRa signal-based outdoor positioning using fingerprint Algorithm[END_REF][START_REF] Carvalho | Hybrid indoor and outdoor localization for elderly care applications with LoRaWAN[END_REF][START_REF] Kim | Feasibility of LoRa for smart home indoor localization[END_REF].

LoRaWAN is a communication protocol and network architecture built on top of the LoRa physical layer. LoRa [START_REF] Vangelista | Frequency shift chirp modulation: The LoRa modulation[END_REF] is a radio technology based on Chirp Spread Spectrum modulation (CSS) operating at industrial, scientific, and medical (ISM) bands (169 MHz, 433 MHz, 868 MHz, 915 MHz, and 2.4 GHz). The most important benefit of the LoRa technology is the CSS modulation which enables a receiver device to demodulate a transmission up to -20 dB below the noise level. Demodulating a LoRa signal below the noise level reduces the battery lifetime of the LoRa device, makes the transmission distance longer, and secures our LoRa communication. All these benefits make LoRa an ideal candidate for indoor tracking. However, detecting a low signal transmission and maintaining a good location accuracy is always a challenge for most indoor tracking applications.

Software Defined Radio (SDR) is a well-known concept in radio tracking, bringing new possibilities, and allowing us to study existing technologies and detect a weak signal transmission and track it. First, a weak signal can be detected by several methods. Detection can be done using some mathematical tools such as Autocorrelation. Autocorrelation [START_REF]Autocorrelation[END_REF] consists of determining the degree of similarity of the same variables between two successive time intervals. Second, tracking the signal can be done using the Angle of Arrival (AoA) positioning and tracking approach. AoA approach has been used in numerous fields, including tracking people and objects.

In this paper, we present an Autocorrelation approach for detecting and measuring the AoA of a weak LoRa signal using an SDR. The main idea consists of detecting the LoRa frequency peaks. The Autocorrelation technique offers an effective way to improve the SNR by reducing the unwanted noise components by a factor of n. Once the signal is detected, The AoA can be calculated using the method described in our previous work [START_REF] Hussein | SYLOIN: Measuring Angle of Arrival of LoRa signals using Software Defined Radio[END_REF]. To the best of our knowledge, no such system can measure the AoA of a LoRa signal up to -20 dB below the noise level.

The rest of the paper is organized as follows. Section II describes the existing LoRa AoA estimation solutions. Section III explains the Autocorrelation method. Section IV presents the implementation and setup of the Fast Autocorrelation method for tracking a LoRa signal using SDR. In section V, we evaluate the effectiveness of the system based on real experiments. In section VI, we conclude the paper and outline future research.

II. RELATED WORKS

The problem of detecting and measuring the AoA of radio signals in localization and tracking services is a rapidly developing area of science and technology. Developing new AoA methods is always a challenge for the research community. LPWAN technologies such as the Long-Range (LoRa) can provide a localization solution based on the time difference of arrival (TDoA) and the received signal strength (RSS) [START_REF] Anjum | RSSI Fingerprinting-Based Localization Using Machine Learning in LoRa Networks[END_REF][START_REF] Lam | RSSI-Based LoRa Localization Systems for Large-Scale Indoor and Outdoor Environments[END_REF][START_REF] Daramouskas | Localization in LoRa Networks Based on Time Difference of Arrival[END_REF][START_REF] Podevijn | TDoA-based outdoor positioning in a public LoRa network[END_REF]. Recently, some research works are trying to mix the RSS, and TDoA techniques with the angle of arrival (AoA) estimation method using the LoRa technology to provide a good localization solution [START_REF] Baik | Hybrid RSSI-AoA Positioning System with Single Time-Modulated Array Receiver for LoRa IoT[END_REF][START_REF] Aernouts | Combining TDoA and AoA with a particle filter in an outdoor LoRaWAN network[END_REF]. LPWAN technologies, including LoRa, utilize a minimum amount of power and a narrow bandwidth to transmit a signal. An RSS-based LoRa localization solution can't offer a good localization accuracy since lowering the power consumption of the LoRa device might result in a lower quality of the RSS and more multipath effect and channel fading. Time-based LoRa localization solutions such as ToA and TDoA also suffer from poor localization accuracy because they require a larger bandwidth or a high timestamp precision of the received signal at the gateway receiver. Addressing these limitations can be done using the AoA estimation technique [START_REF] Bnilam | Adaptive probabilistic model using angle of arrival estimation for IoT indoor localization[END_REF]. Using an AoAbased localization solution has several advantages. First, unlike the TDoA technique, the AoA technique does not require a large signal bandwidth and a high timestamp precision at the receiver. Secondly, the AoA-based localization solutions measure the angle between the transmitter and the receiver by measuring the phase of the received signal at the receiver gateway. Consequently, the narrow bandwidth of the LoRa signal will not affect the phase of the signal and, therefore, the accuracy of the AoA of the received signal. Third, the AoA technique does not depend on the received signal strength and, therefore, the phase and the angle of the received signal will not be affected by the multipath and channel fading. These benefits make the AoA a promising localization solution for the indoor environment.

Recently, AoA-based LoRa localization solutions have been introduced in the literature. In [START_REF] Bnilam | LoRay: AoA estimation system for long range communication networks[END_REF], an AoA localization technique of LoRa signals based on an SDR receiver has been proposed. The proposed method found an estimation error of 2 • in LOS (Line-of-Sight), and 10 • in NLOS. However, this method consists of 8-element uniform linear array (ULA). The 8 ULA consists of 8 independent SDR receivers for which a hand-crafted crystal oscillator is needed to achieve synchronization between the channels. Another AoA-based LoRa localization approach based on phase detectors has been presented in [START_REF] Bnilam | Angle of Arrival Estimation System for LoRa Technology based on Phase Detectors[END_REF]. The phase detectors make a cost-effective solution to measure the AoA with median and maximum estimation errors of 4 • and 14 • , respectively.

Recently, we introduced the SDR-based localization system that can estimate the AoA of a LoRa signal by measuring the phase difference of the received signal at two antennas [START_REF] Hussein | SYLOIN: Measuring Angle of Arrival of LoRa signals using Software Defined Radio[END_REF]. We built our system based on an individual universal software radio peripheral (USRP) B210 SDR receiver. The choice of this front-end receiver was motivated by its two coherent RX channels. The prototype has been tested in an indoor laboratory office environment a with lot of furniture. The experimental results show that the AoA of the LoRa signal can be measured accurately with a maximum estimation error of 5 • in the range of [0 • ,180 • ]. However, our system can only achieve an accurate AoA measurement with high signalto-noise ratio (SNR) communication. In this paper, we aim to validate our AoA-localization system with a very low SNR up to -20 dB. A new approach is proposed to detect and measure the AoA of a very weak LoRa transmission signal. The presented approach is based on two-antenna reception to exploit the Fast autocorrelation signal detection and the phase difference AoA measurement. We, therefore, believe that an accurate AoA of a weak LoRa signal reduces the battery lifetime, extends the communication lifetime, and secures the communication of the LoRa devices. To the best of our knowledge, we are the first to implement a real AoA for a very weak LoRa signal.

In the following section, we will provide an overview of the two-antenna autocorrelation technique for weak signal detection.

III. DUAL-ANTENNA FAST-AUTOCORRELATION

A. Autocorrelation Function

An Autocorrelation is a signal processing tool used to separate useful signal from noise. The Autocorrelation analysis consists of finding the presence of useful segments in a signal. The autocorrelation of a signal is just the correlation coefficient. However, instead of applying the correlation between two different signals, the correlation is between the two segments of the same signal at times T i and T i+k .

Suppose we have two continuous signals x(n), y(n), n = 1,2,3,..,N, where N is the number of samples. The discrete correlation function of the two signals can be calculated by the following equation:

Corr(x, y) k = 1 N N -1 n=0 x k+n y n (1)
Equation ( 1) measures the correlation function of the two signals. When the two signals are similar, then their correlation coefficient will be 1. However, if the two signals are different, then their correlation coefficient will be 0. If we suppose that = y(n) = x(n), then Equation (1) will be:

Corr(x, x) k = 1 N N -1 n=0 x k+n x n (2)
Equation ( 2) measures the correlation function between two same signals x(n) and x(n+k) separated by a time lag k, which is the autocorrelation function of the signal x(n). The autocorrelation function presented by (2) reaches its maximum value at the lag=0, and is always symmetric around zero, which means that the autocorrelation values are equals at lags +k and -k. Using these properties, it is possible to detect weak signals. Using the FFT with the autocorrelation function is described in the following section.

B. Fast Fourier Transform Autocorrelation

An FFT is an important signal processing tool used in most digital signal processing (DSP) applications. As the name implies, the FFT is a fast version of the Discrete Fourier Transform (DFT). FFT is an algorithm used to perform a series of computations that computes DFT much more rapidly than the other available algorithms. Therefore, using FFT will reduce the number of calculations necessary. Let x(n) be a discrete signal, the Discrete Fourier Transform of x(n) is given by:

X k = N -1 n=0 x n e -i2πkn N k = 0, 1, 2, .., N -1 (3)
Equation ( 3) represents the DFT computation of N-point data sequence with O(N 2 ) multiplies and adds. Thus, when using FFT, the computation would take O(N logN ) multiplies and adds.

Signal detection is a key function of the FFT. Using FFT and the autocorrelation function can accelerate the autocorrelation computation. The Wiener-Khinchin theorem [START_REF] Bnilam | Angle of Arrival Estimation System for LoRa Technology based on Phase Detectors[END_REF] states that the autocorrelation function and the Power Spectral Density (PSD) of the signal are Fourier pairs. The autocorrelation function is the Fourier transform of the PSD.

Let x(n), y(n) be two discrete signals, the correlation function is defined by:

Corr(x, y) k ⇔ X(f )Y * (f ) (4)
From Equation (4), the correlation function of the two signals can be obtained by multiplying the FFT of one signal by the complex conjugate of the FFT of the other signal. As the autocorrelation function of x(n) is the correlation with its delay copy, the autocorrelation will therefore be as follows:

Corr(x, x) k ⇔ X(f )X * (f ) ≡ |X(f )| 2 (5)
Where |X(f )| 2 is the Power Spectral Density of the signal. Using these properties of autocorrelation fucntion and FFT, it is possible to detect weak signal buried in noise.

The next section describes how to detect weak signals using Autocorrelation, FFT, and Dual-Antenna detection.

C. Dual-Antenna weak signal detection

This section discusses the Dual-Antenna Autocorrelation detection. The proposed detection method measures the Autocorrelation function between outputs of two sensors of the receiving antenna array.

The dual-Antenna approach is a very promising topic in radio communication systems. For the receiver side, Dual-Antenna offers a high potential solution to receive two replicas of the primary signal. The Correlation method can be used to determine whether the primary signal exists or not.

Figure 1 shows two signals X1(t) and X2(t) received by the antennas A1 and A2 respectively. When the transmitter and the receiver are in Line-Of-Sight (LOS) conditions, there will be a time delay ∆(t) between the two signals received at the two antennas. Since the two signals X1(t) and X2(t) are the same time series data received at a different time by the two antennas, X1(t) = X2(t + ∆(t)).

Using the properties of the correlation function between two same signals separated by a time lag t given by the Equation (2), and the Fast Autocorrelation function using FFT given by Equation ( 5), the detection of the signal using a Dual-Antenna receiver is possible.

By collecting enough samples of the received signals, the Fast Autocorrelation, of the signals received from the two Antennas can be built. Based on the Autocorrelation peaks finding and the threshold value, signal detection can be carried out.

In the next section, we present our solution of detecting a LoRa signal using a Dual-Antenna USRP B210 SDR Receiver with the help of GNU Radio.

IV. PROPOSED METHOD

The following method consists of detecting weak LoRa signals on a Dual-Antenna Receiver by measuring the Fast Autocorrelation function. We first determine whether the LoRa signal exists by using the Fast Autocorrelation signal detection. Then, we measure the AoA of the detected signal based on our published work [START_REF] Hussein | SYLOIN: Measuring Angle of Arrival of LoRa signals using Software Defined Radio[END_REF]. Finally, we compare the two methods in terms of accuracy level. The USRP B210 SDR receiver and the GNU Radio development environment form the basis for this research work.

A. USRP-Based SDR Receiver

Most Traditional radios use integrated circuits and specialized equipment for transmission and reception. Traditional radios are made up of mixers, amplifiers, analog/digital converters, modulators, and demodulators. An innovative concept named software defined radio (SDR) replaced the analogical traditional radios. In an SDR, signal processing functions are performed by a single device, which performs the functions of multiple hardware radios. Most of the functionality of the SDR are implemented in software rather than in hardware. GNU Radio founder Eric Blossom cites that the SDR is the technique of getting code as close as possible to the antenna. It turns radio hardware problems into software problems [START_REF] Eric | Exploring Gnu Radio[END_REF].

Ettus Research is one of the companies that have developed an open source hardware platform dedicated to the SDR. The development was focused on the Universal Software Radio Peripheral (USRP). The USRP is characterized by its accuracy, large bandwidth and low price compared to other hardware systems. This open source device defines the interface between the RF analog domain and the General Purpose Processor (GPP) and can generate and receive any kind of signals. The USRP can receive radio waves captured by the antenna and convert these waves into a digital signal that can be processed by the computer using an analog to digital converter (ADC). On the other hand, to generate a radio wave, the USRP converts a digital wave processed by the computer into a RF signal using a digital to analog converter (DAC).

The model of USRP used in this research was the USRP B210. The USRP B210 is a Dual channel SDR receiver which contains 2TX and 2RX ports with fully coherent 2 x 2 MIMO (Multiple Input Multiple Output) capabilities. The coherency between the 2RX ports is an important feature for computing the phase difference and the AoA between the received signals [START_REF] Hussein | SYLOIN: Measuring Angle of Arrival of LoRa signals using Software Defined Radio[END_REF]. The B210 covers frequencies from 70 MHz to 6 GHz with a maximum bandwidth of 56 MHz in real-time. Moreover, this SDR board combines the Spartan6 FPGA, and the AD9361 transceiver that performs the conversion of ADC and DSP. The B210 can be connected directly to a computer via a high-speed Ethernet local area network cable or a USB 2.0/3.0 cable and is compatible with the software development platform GNU Radio.

B. GNU Radio

GNU Radio is one of the most popular software platforms used to interface with a USRP SDR. The USRP implements the hardware part of the system, whereas GNU Radio handles the software part.

GNU Radio provides signal processing blocks for the development of SDR applications. The blocks are written using the C++ programming language, and the connection between them is provided by the Python programming language. GNU Radio is widely used by academic research as well as in the industrial and commercial domain.

In the rest of this paper, the USRP B210 is used combined with GNU Radio development software platform.

C. USRP B210 Fast Autocorrelation signal detection

The block diagram in Figure 1 shows the basic principle of our proposed solution. The process starts by extracting the two signals received at the two channels using a UHD: USRP Source block which receives the data from the B210 receiver through the USB3.0 cable. Thereafter, the two signals diverge into two paths. For each path, the output of the UHD: USRP Source block is passed through a Stream to Vector block. The Stream to Vector block transforms the data in I and Q format. Now that the data are in I and Q format, the FFT block computes the FFT to transform the data from time domain to frequency domain. Up until this point the data are in the Once the maximum of the peak value is found, the output of the Max block is passed to the Threshold block. The Threshold block is used to compare the output of the Fast Autocorrelation detection algorithm with a predefined threshold. This block has two parameters: High value and Low value. If the output of the Max block is higher than the value set for the High parameter, the output of the The Threshold block is 1. On the other hand, if the output of the Max block goes below the value set for the Low parameter, the output of the Threshold block is 0. The decision output of the Fast Autocorrelation detection algorithm is given by the following equation:

Decision = 1 if M ax > α 0 if M ax ≤ α (6)
where M ax is the peak of the Fast Autocorrelation detection and α is the predefined threshold. The predefined threshold α is determined using a post processing analysis. The F ileSink block is used to write the data in an empty file. The written data are processed in Matlab to determine the optimum predefined threshold value.

The case when the value of the Decision is equal to 1 meaning that the signal is detected, and the case when the value of the Decision is equal to 0 meaning that there is no signal. Once the signal is detected, the algorithm computes the AoA using the phase difference method implemented in [START_REF] Hussein | SYLOIN: Measuring Angle of Arrival of LoRa signals using Software Defined Radio[END_REF]. The implementation of each block shown in Figure 1 is shown in Figure 2. A full description of the phase difference and AoA measurement blocks can be found in [START_REF] Hussein | SYLOIN: Measuring Angle of Arrival of LoRa signals using Software Defined Radio[END_REF].

V. EXPERIMENTATION AND RESULTS

The proposed method was evaluated by performing an experimental measurement in our indoor research laboratory. The equipment setup and experimental measurement results will be discussed in detail in the following sections.

A. Equipment Setup

The FFT Autocorrelation signal detection and AoA measurement is implemented with GNU Radio and USRP B210 SDR receiver for capturing a LoRa signal.

The real time LoRa signal captured by the USRP B210 is transmitted by the Adafruit Feather M0 transceiver containing a RFM95 LoRa module. The Adafruit Feather M0 combines an Arduino-based microcontroller with a LoRa radio module. The radio transceiver can transmit or receive radio signals at both frequencies 868 MHz and 915 MHz, which can be programmed in software. The Feather M0 transceiver was programmed using the Arduino IDE development toolkit via a USB-to-Serial program. The transmitted signal was configured with the frequency 868 MHz.

In the receiving end, the signal was received using the USRP B210 SDR which was connected using USB3.0 high speed connectivity to a PC running Linux. The SDR card was configured with the GNU Radio Companion (GRC) 3.7 and the USRP hardware Driver (UHD) 3.14.0. The GNU Radio Companion (GRC) is a graphical user interface used to compile GNU Radio flow graphs. The B210 was configured with the following parameters: center frequency: 868 MHz, sampling rate = 2 MHz, and the bandwidth is 1 MHz.

The Feather M0 LoRa transceiver and the USRP B210 are both equipped with a quarter-wavelength commercial off-theshelf omnidirectional antenna. Both transmitter and receiver were fixed on a tripod with a Line-Of-Sight (LOS) condition. The complete equipment setup is shown in Figure 3.

B. Measurement data

For the measurement data, a LoRa signal was transmitted using a Feather M0 with the following parameters: frequency = 868 MHz, bandwidth = 125 kHz, spreading factor = 7, and coding rate = 4/5. On the USRP B210 side, we assume that the sampling frequency is equal to 2 MHz, and the number of sampled points (FFT length) is equal to 1024 by default.

The experiments were divided into two parts: the first was to identify the suitable FFT length in order to improve the detection accuracy at very low power transmission (TX) or signal to noise ratio (SNR). The second was to measure the AoA of the LoRa signal at each SNR level.

In experiment 1, the transmitter is programmed to transmit a signal with a transmission power going from +20 dBm to -1 dBm. The detection accuracy of the proposed method is measured for different FFT length. The performance of the algorithm is shown in Figure 4. When the FFT length is equal to 1024, the LoRa signal can be detected effectively when the SNR is greater than -12 dB. When the FFT length is equal to 65536, the detection of the LoRa signal increases to -20 dB. The detection performance of the LoRa signal can be improved effectively when increasing the FFT length because the autocorrelation of noise become smaller when increasing the FFT length. Figure 5 and6 show respectively the values of the autocorrelation peaks for the LoRa signal with two different FFT length. It clearly appears that when using an FFT length equal to 65536 (sample size 30 Hz), the autocorrelation peaks are easier to detect than using an FFT length equal to 1024 (sample size 2 kHz). Therefore, the detection performance of the algorithm is better at high FFT length. In our next experiment, the FFT length is chosen to be equal to 65536.

In experiment 2, we measured the AoA of the LoRa signal at three different SNR levels: 20 dB, 0 dB, and -20 dB. For each SNR level, the AoA was measured in the range of [0, 180] degrees with 15 degrees increments. The results of this measurement is illustrated in Figure 7. We can clearly notice that for an SNR equal to 20dB, the measured AoA are very close to the real AoA values. The angular deviation increases significantly when the measurements were taken with SNR equal to 0dB. For an SNR equal to -20dB, the algorithm still gives a good measurements accuracy for AoA between 45 degrees and 135 degrees. The angular deviation increases exponentially for AoA less than 45 degrees and greater than 135 degrees. Therefore, the AoA accuracy is better for large angles when the SNR is equal to -20 dB. The AoA deviation error of this experiment is shown in Figure 8. When the SNR is equal to 20 dB, the algorithm can successfully measure the AoA with a deviation error less than 4 degrees in the range of [45,135] degrees, and less than 6 degrees for small angles. Otherwise, when the SNR decreases to 0 dB and -20 dB, the deviation error increases to 9 and 12 degrees respectively for small angles.

According to these experimental results, we can indicate an efficient method to deploy an indoor location infrastructure exploiting low AoA. Indeed, the alignment of the dual antennas of a LoRa gateway must favor the angles located between 45 and 135 degrees in the building so that the other angle values are as much as possible outside the building as can be seen on the Figure 9. Thus, the accuracy of the AoA will be in the range 7 to 8 degrees. With this topology, the best accuracy will be obtained at the center of the building where the AoA are in the green area, while the worst accuracy will be obtained at the building edges where the AoA are in the red area. For an angle with 8 degrees of accuracy, the accuracy will be about 1 meter at a distance of 8 meters. To maintain an accuracy of less than one meter, the distance between two LoRa gateways in the building will be approximately 16 meters. During a critic situation, it will be possible for a gateway to send a specific message to an object in order to ask for a full transmit power reply to improve the accuracy with a factor of nearly 4. An accurate AoA on weak signals allow us to strengthen our existing localization algorithms which will merge signal attenuation, time of flight and AoA to make our coordinate calculations more robust. Moreover, our location functions are able to work in a larger space. Additionally, the asymmetry of the transmission power between the gateway and the object makes it possible to extend the communication lifetime of the objects both by reducing the transmission power and by facilitating the reception of strong signals coming from the gateway. Allowing objects to transmit weakly improves the stealth and the confidentiality of exchanges and makes it more difficult for an enemy device to receive and locate our broadcasting. Using low power transmissions allow us a significant budget link inside a building because the reception threshold is very low, between -125 and -140 dBm, conducting to a communication distance longer than Wi-Fi with a power transmission at 100 mW. Finally, we believe that it will be possible to decode the simultaneous transmission of 2 or 3 objects to avoid collisions and limit re-transmissions which would become unnecessary.

Fig. 1 .

 1 Fig. 1. A GNU Radio flow graph of the proposed system

Fig. 2 .

 2 Fig. 2. A GRC flow graph for the USRP B210 to compute the FFT Autocorrelation and AoA Measurements

Fig. 3 .

 3 Fig. 3. The complete installation setup. Both Feather M0 and USRP B210 are fixed on the tripod.

Fig. 4 .Fig. 5 .Fig. 6 .

 456 Fig. 4. Autocorrelation detection accuracy for different FFT length

Fig. 7 .Fig. 8 .Fig. 9 .

 789 Fig. 7. AoA measurements with different SNR level

ACKNOWLEDGMENT

This work is a partnership with Jolliet and SECA companies and is funded by the Franche-Comté region, the European FEDER fund and the French Public Investment Bank (Banque Publique d'Investissement). Authors want to warmly thank Stephane Givron for his highly valuable contribution.