A Simon

R Pascaud

T Callegari

L Liard

O Pascal

O Pigaglio

Static and Dynamic Control of Limiting Threshold in Plasma-based Microstrip Microwave Power Limiter

Keywords: microwave limiters, radio frequency (RF) tunability, plasma devices, gas discharge devices

Static and dynamic control of the limiting threshold in a very low-loss plasma-based microstrip power limiter is investigated in order to prevent receivers from being threatened by high-power microwave (HPM). An analytic model of the microstrip circuit is proposed to derive the influence of its design on the limiting threshold of the microwave power limiter. Experimental results are in good agreement with those expected. Finally, an original approach to allow discrete limiting threshold tunability is presented and validated experimentally.

I. INTRODUCTION

Current receiver systems can be highly vulnerable to highpower microwave (HPM) [START_REF] Backstrom | Susceptibility of electronic systems to high-power microwaves: summary of test experience[END_REF]. To protect them from these threats, microwave power limiters are commonly involved. They usually rely on PIN diodes, MEMS, and other non-linear devices [START_REF] Yang | A novel analysis of a ku -band planar p-i-n diode limiter[END_REF], [START_REF] Shojaei-Asanjan | Tunable rf memsbased frequency dependent power limiter[END_REF]. However, when the threat becomes important these components deteriorate and lose their power limiting properties. Recently, limiters based on plasma discharges have been proposed to handle very high-power threats.

In microwave power limiters, plasma discharges arise when the incoming HPM is sufficient to ionize an inert gas. The nonlinear interaction between the incident electromagnetic wave and the ionized gas is finally used for power limitation. Plasma discharges have been widely used in T/R (Transmit/Receive) tubes in order to protect RADAR receivers [START_REF] Smullin | Microwave Duplexers[END_REF]. If T/R tubes are fully compatible with waveguide technology, they can not be easily integrated into compact microstrip circuits.

Several works have thus studied the integration of plasma discharges into planar microwave power limiters. The first proposed solution was to integrate a plasma discharge into a microstrip line [START_REF] Patel | Microstrip plasma limiter[END_REF]. The measured limiting threshold of this device was above 50 W which remains important. In order to reduce the limiting threshold, recent works have considered resonant circuits [START_REF] Cross | Theory and demonstration of narrowband bent hairpin filters integrated with accoupled plasma limiter elements[END_REF], pre-ionization techniques [START_REF] Pascaud | Plasma microdischarge as power-induced limiter element in microstrip devices[END_REF], or both [START_REF] Semnani | Power limiting characteristic of a plasma-loaded evanescent-mode cavity resonator[END_REF], [START_REF] Pascaud | Low insertion loss microplasma-based limiter integrated into a microstrip bandpass filter[END_REF]. Using such techniques, the limiting threshold has been reduced from tens of watts to hundreds of milliwatts. However, no design guidelines have been provided to tune the limiting threshold.

In this paper, we study both theoretically and experimentally the control of the limiting threshold of a plasma-based microwave power limiter by using several geometrical parameters of the microstrip circuit. Besides, an original, simple, and low-loss technique is presented to provide tunable limiting threshold.

II. BASIC PRINCIPLE

In plasma-based power limiters, the limiting capability is achieved when the gas breakdown happens. It practically corresponds to a critical value of incident microwave electric field |E B | that only depends on plasma parameters related to ionization processes such as the type of gas, the electrode material, or pre-ionization. However, from a circuit point of view, the electric field magnitude at any point is directly linked to the incident microwave power, but also to the microwave circuit properties.

In order to illustrate the influence of the microwave circuit design on the limiting threshold, we here consider a microstrip annular ring resonator. Nevertheless, the proposed analysis remains valid for any microwave circuit. According to [START_REF] Hsieh | Equivalent lumped elements G, L, C, and unloaded Q's of closed-and open-loop ring resonators[END_REF], the input impedance of a low-loss unloaded microstrip ring resonator at the resonant frequency of its first mode is:

Z in Z 0 Q unloaded 2π (1)
Considering the cos(ϕ) angular dependence of the electric field for the fundamental mode, and since the voltage at the feeding point is V in = √ 2P in Z in , we can show that the magnitude of the electric field along the ring resonator may be written as:

|E| = 1 h ms Z 0 Q unloaded π P in | cos(ϕ)| (2)
with h ms the substrate thickness, Z 0 the characteristic impedance of the ring resonator, Q unloaded its unloaded quality factor, P in the input power, and ϕ the angle around the ring. Equation [START_REF] Yang | A novel analysis of a ku -band planar p-i-n diode limiter[END_REF] shows that modifying h ms , Z 0 , or ϕ can theoretically change the required input power to reach |E B | and thus control the limiting threshold of a plasma-based microwave power limiter.

III. DEVICES AND EXPERIMENTAL SETUP In order to experimentally verify the control of the limiting threshold, several plasma-based microwave power limiters have been designed according to the schematic presented in Fig. 1a. It consists of an electromagnetically coupled annular ring resonator in microstrip technology combined with a plasma pre-ionization system [START_REF] Pascaud | Low insertion loss microplasma-based limiter integrated into a microstrip bandpass filter[END_REF]. The pre-ionization is ensured by a micro hollow cathode discharge (MHCD) that has been optimized to work at a pressure of 10 Torr (i.e., 1333 Pa) in argon that is to say h MHCD = d MHCD = 1.5 mm.

Here, the proposed annular ring resonator is loaded at its input and output. This is obviously different from the case considered in Section II but if the loading and the coupling modify the electric field magnitude in the ring, it can still be assumed that this electric field has the same dependency on h ms , Z 0 , and ϕ as stated in (2). A dedicated experimental setup has been developed to perform a power budget for different configurations. A microwave unit allows to set up a power ramp and measure the input and output power at the resonant frequency of each resonator. Moreover, the experimental setup includes additional units that control the type of gas, its pressure, and the DC pre-ionization.

IV. STATIC CONTROL OF LIMITING THRESHOLD BY CIRCUIT DESIGN

A. Circuits design

In order to illustrate the static control of the limiting threshold, three ring resonators have been designed using fullwave simulations. Each ring resonator has been realized on a Rogers RO4003C substrate (ε r = 3.55) and optimized to have a resonant frequency equal to 2.5 GHz. A reference ring resonator (Resonator 1) has been printed on a 0.508 mm thick substrate with a 50 Ω characteristic impedance for the ring resonator and a hole located at ϕ d = 0 • (see Fig. 1a), where the magnitude of the electric field reaches a maximum at resonance. Then two additional ring resonators with each a different parameter compared to Resonator 1 have been designed. The first one has a 25 Ω characteristic impedance for the ring (Resonator 2), while the other one has a 1.524 mm substrate thickness (Resonator 3). Table I summarizes the dimensions of each microstrip circuit. Note that in each case the coupled microstrip lines have a 50 Ω characteristic impedance.

Fig. 2 presents the simulated and measured S-parameters for each resonator. One can distinguish a slight frequency shift between numerical simulations and measurements mainly due to the manufacturing tolerances. These simulations were performed with ANSYS HFSS.

B. Experimental results and discussion

Fig. 3 presents the output power versus the input power for each resonator at the resonant frequency with or without preionization. Measurements using a VNA have shown that the additional insertion loss due to pre-ionization remains below As shown in Fig. 3, without pre-ionization, the three resonators remain in a linear regime. This microwave breakdown happens for an input power larger than 40 dBm. With preionization, the limiting threshold values are well lowered, and the typical response of a power limiter is observed with a linear regime at low input power and a saturated behaviour at high input power.

The limiting threshold is defined as the 1 dB compression point. Resonators 1, 2, and 3 have a measured limiting threshold of 16.2 dBm, 19.6 dBm, and 24.8 dBm, respectively. Considering (2), and assuming that the required breakdown electric field is constant in each case, the theoretical differences between the limiting thresholds taking Resonator 1 as a reference are 3 dB for Resonator 2, and 9.5 dB for Resonator 3. The comparison of the relative limiting thresholds with the experimental results shows an accuracy of 1 dB. This difference is mainly attributed to the measurements accuracy and the different assumptions that have been made to establish the theoretical model. However, these results clearly exhibit the capability of having a static control of the limiting threshold for a given power limiter.

V. DYNAMIC CONTROL OF LIMITING THRESHOLD WITH MULTIPLE PLASMA DISCHARGES

A. Circuit design

If h ms and Z 0 are fixed once the circuit has been realized, one can easily use multiple holes in the microstrip circuit in order to dynamically control the position of the discharge ϕ d . Such a property paves the way of discrete limiting threshold tunability. Thus, an additional ring resonator including three holes drilled at ϕ d1 = 0 • , ϕ d2 = 30 • , and ϕ d3 = 60 • has been printed on a 0.508 mm thick Rogers RO4003C substrate. The dimensions of this resonator are similar to those of the Resonator 1 (see Table I). In order to control independently the pre-ionization of each plasma hole, a new MHCD circuit was printed on a 1.524 mm thick Rogers RO4003C including three cathodes at positions ϕ d1 , ϕ d2 , and ϕ d3 .

B. Experimental results and discussion

This microstrip ring resonator has been characterized at its resonant frequency by pre-ionizing independently each MHCD. Fig. 4 presents the output power versus the input power at the resonant frequency with or without pre-ionization. Again, a current of 0.5 mA is considered here. In this case, we can observe that the transition between the linear and the non-linear behaviour is abrupt for the ϕ d2 = 30 • and ϕ d3 = 60 • positions contrary to the smooth transitions seen in Fig. 3. This can be due to the observed shift of the plasma discharge from the ϕ d2,3 positions (Picture 1) to the ϕ d1 position (Picture 2). It can be explained by the fact that the microwave plasma discharge, once ignited, goes where the microwave electric field is maximum. Regarding the ϕ d1 case, it should not however have this abrupt behaviour. We suppose that there were some issues with the MHCD and microstrip holes alignment within the manufacturing process.

Fig. 4 shows that the limiter turns on when the input power exceeds 20.5 dBm for ϕ d1 , 22.1 dBm for ϕ d2 , and 26.3 dBm for ϕ d3 . Considering (2), and assuming again that the required breakdown electric field is constant in each case, the theoretical differences between the limiting thresholds taking the ϕ d1 case as a reference are 1.2 dB for ϕ d2 and 6 dB for ϕ d3 . The measured thresholds show a difference compliant with the theoretical results and prove that this method can be considered as accurate to enable dynamic control of the limiting threshold by pre-ionizing the appropriate hole.

VI. CONCLUSION

Static and dynamic control for microstrip power limiters have been theoretically and experimentally investigated. The theoretical method used in this paper exhibits a good accuracy that can be employed to design any microstrip plasmabased power limiter. Besides, an original approach to obtain discrete limiting threshold tunability has been proposed and experimentally validated. Such a low-loss solution may be very promising to tune the level of receiver protection depending on the electromagnetic environment.

VII. ACKNOWLEDGEMENT

Fig. 1 .

 1 Fig. 1. Integration of a plasma discharge in a microstrip circuit

Fig. 2 .

 2 Fig. 2. Simulated and measured S parameters of each resonator

Fig. 3 .

 3 Fig. 3. Measured output power versus input power for the three microstrip ring resonators

Fig. 4 .

 4 Fig. 4. Measured output power versus input power of the microstrip resonator with different plasma discharge positions.

The authors would like to thank the support from the Université Fédérale Toulouse-Midi-Pyrénées (UFTMiP) and the Région Occitanie. This work was also supported in part by the French National Research Agency under Grant ANR-15-ASTR-0020-02, DIOMEDE project, within the framework of the 2015 ASTRID program.


Daniel Quemada
email: danielquemada@orange.fr

Claudio Berli

DESCRIBING THE FLOW CURVE OF SHEAR-BANDING FLUIDS THROUGH A STRUCTURAL MINIMAL MODEL

Keywords: Shear-banding, Couette rheometry, wormlike micellar solutions, structural model

Main characteristics of colloidal systems that develop fluid phases with different mechanical properties, namely shear-banding fluids, are briefly reviewed both from experimental and theoretical (modelling) point of view. A non-monotonic shear stress vs. shear rate constitutive relation is presented. This relation derives from a phenomenological model of a shear ratedependent viscosity describing structural changes and involves the possibility of multivalued shear rates under a given shear stress. In the case of a stress-dependent viscosity, the same model allows one to predict vorticity banding. Predictions of this model under controlled stress are discussed, namely occurrence of a kind of top-and bottom-jumping of the shear rate in response to stress increasing-decreasing. Applying this model to evaluation of the flow curve of such colloidal systems is performed. Particular emphasis is placed on the adequate computation of the shear rate function in cylindrical Couette cells in order to handle the corresponding flow curve which exhibits the well-known shear stress plateau. Indeed, as different fluid phases coexist in the flow domain, measured (torque vs. angular velocity) data cannot be directly converted into rheometric (shear stress vs. shear rate) functions. As the lacking non-local terms in the model prevents the direct determination of the stress-plateau, this value is included as an adjustable parameter. Thus model predictions satisfactorily match up experimental data of wormlike micellar solutions from the literature.

-Introduction

Under shear flow conditions, many complex fluids develop fluid phases with different mechanical properties, namely shear bands, parallel to the flow direction, where the shear rate takes different mean values. This so-called shear banding behaviour has been interpreted as resulting from flow instabilities associated to the existence of a non-monotonic flow curve with a region of negative slope. Under controlled stress experiments, band coexistence occurs in a certain range of shear rate where the flow curve presents either a stress plateau or a notably reduced slope. Typical results found in the analysis of wormlike micellar solutions (WMS) are represented in Fig. 1.

Shear banding has mainly been found experimentally in WMS [START_REF] Berret | Isotropic-to-nematic transition in wormlike micelles under shear[END_REF][START_REF] Grand | Slow Transients and Metastability in Wormlike Micelle Rheology[END_REF][START_REF] Britton | Two-Phase Shear Band Structures at Uniform Stress[END_REF][START_REF] Lerouge | Correlations between Rheological and Optical Properties of a Micellar Solution under Shear Banding Flow[END_REF][START_REF] Berret | Rheology of Wormlike Micelles: Equilibrium Properties and Shear Banding Transition[END_REF][START_REF] Cappelaere | Rheology, Birefringence, and Small-Angle Neutron Scattering in a Charged Micellar System: Evidence of a Shear-Induced Phase Transition[END_REF], colloidal crystals [START_REF] Chen | Rheological and Microstructural Transitions in Colloidal Crystals[END_REF][START_REF] Cohen | Slip, Yield, and Bands in Colloidal Crystals under Oscillatory Shear[END_REF][START_REF] Bonn | Bistability in non-Newtonian flow: Rheology of lyotropic liquid crystals[END_REF] and lamellar surfactant systems [START_REF] Diat | Effect of shear on a lyotropic lamellar phase[END_REF][START_REF] Salmon | Shear banding in a lyotropic lamellar phase. I. Timeaveraged velocity profiles[END_REF]. In these materials, shear-induced changes of the microstructure have been observed by birefringence [START_REF] Lerouge | Correlations between Rheological and Optical Properties of a Micellar Solution under Shear Banding Flow[END_REF][START_REF] Cappelaere | Rheology, Birefringence, and Small-Angle Neutron Scattering in a Charged Micellar System: Evidence of a Shear-Induced Phase Transition[END_REF][START_REF] Salmon | Shear banding in a lyotropic lamellar phase. I. Timeaveraged velocity profiles[END_REF], local light and neutron scattering [START_REF] Makhloufi | Rheo-optical Study of Worm-like Micelles Undergoing a Shear Banding Flow[END_REF][START_REF] Dhont | Shear-banding and microstructure of colloids in shear flow[END_REF], associated with various methods of velocimetry [START_REF] Salmon | Velocity Profiles in Shear-Banding Wormlike Micelles[END_REF][START_REF] Salmon | Shear banding in a lyotropic lamellar phase. I. Timeaveraged velocity profiles[END_REF][START_REF] Britton | Shear banding instability in wormlike micellar solutions[END_REF][START_REF] Fischer | Shear banding and the isotropic-to-nematic transition in wormlike micelles[END_REF][START_REF] Méndez-Sánchez | Particle image velocimetry of the unstable capillary flow of a micellar solution[END_REF][START_REF] Miller | Transient evolution of shear-banding wormlike micellar solutions[END_REF], performed in various geometries (pipe, cone-and-plate and coaxial cylinders). In many cases, experimental findings prove that banding can be interpreted as a I-N phase transition between a disordered (Isotropic) phase coexisting with an ordered (Nematic) one. Hence the fluid exhibits shearthinning, with low viscosity as micelles align in the flow direction at high shear rate. By contrast, micelle entanglements lead to high viscosity at low shear rate. In Fig. 1, inset boxes schematically show the microstructural changes driven by shear flow. General agreement between birefringence and velocity profiles have been observed, although some systems exhibit a high viscous nematic state at high shear rate, which could reveal either a possible mesoscale ordering [START_REF] Fischer | Shear banding and the isotropic-to-nematic transition in wormlike micelles[END_REF] or shear-thickening effects as yet observed in similar systems in this rate domain [START_REF] Hartmann | Shear Thickening of an Aqueous Micellar Solution of CTAB and Sodium Tosylate[END_REF].

Another kind of banding formation is called vorticity banding because band interfaces appear as turbid and clear (pancake-like) rings oriented perpendicular to the vorticity axis. However no steady state has been observed in this case due to stress oscillations that take place once the pancake structure appears [START_REF] Fischer | Time dependent flow in equimolar micellar solutions: transient behaviour of the shear stress and first normal stress difference in shear induced structures coupled with flow instabilities[END_REF][START_REF] Fischer | Shear-banding structure orientated in the vorticity direction observed for equimolar micellar solution[END_REF][START_REF] Herle | Alternating Vorticity Bands in a Solution of Wormlike Micelles[END_REF].

A large number of theoretical works have demonstrated that the main feature in shear banding modelling is the need of a non-monotonic constitutive relation in order to obtain a flow curve, shear stress σ vs. shear rate γ& , which exhibits a region where the slope

γ σ & d d is negative.
Basically, it is admitted that the flow curve should take the form shown in Fig. 2, with maximum and minimum shear stresses, σ M and σ m , reached at shear rates M γ& and m γ& , respectively. For any given stress σ 0 lying within the interval [σ m , σ M ], one may define two "stable" shear rates, ) requires that σ remains uniform across the gap, σ 0 = σ( 1 γ& ) = σ( 2 γ& ) , which suggests that the flow is shared out in two coexisting bands with constant rates 1 γ& and 2 γ& . The average shear rate thus follows a lever rule

2 1 av ) x 1 ( x γ γ γ & & & - + =
, where x gives the interface location within the gap of the rheometric cell, with 0 < x < 1 [START_REF]these factors x and (1-x) being sometimes identified to the volume fractions of each band[END_REF]. All streamlines being equivalent in steady planar Couette flow, the presence of shear bands should be possible for any stress σ 0 such as σ m <σ 0 < σ M , thus band interface location cannot be determined. However, experiments in real systems have shown the robustness of a well defined stress plateau, σ = σ *, at which the two phases coexist. Some criterion is required for selecting this value.

Therefore, modelling should incorporate some new ingredient. Note that in the case of nonplanar geometries, the stress variation within the gap cannot give the selection criterion, although it might explain the experimental observation of a quasi-plateau in this curve instead of a true one.

Although the way of spontaneous formation of a banded flow still remains an open question, several works have proposed that instabilities should result from a flow/structure coupling, i.e. coupling between shear induced changes of the internal structure of the material (as the isotropic-to-nematic transition in WMS) and pure hydrodynamic instabilities associated to the presence of bands which induce non-homogeneous velocity profile [START_REF] Berret | Isotropic-to-nematic transition in wormlike micelles under shear[END_REF][START_REF] Fischer | Time dependent flow in equimolar micellar solutions: transient behaviour of the shear stress and first normal stress difference in shear induced structures coupled with flow instabilities[END_REF].

The basic role played by molecular orientation in I-N transition and more generally in nematic structures has been treated by Doi and Hess theories which appear as the more powerfull tools to solve this problem [START_REF] Hess | Fokker-Planck-equation approach to flow alignment in liquid crystals[END_REF][START_REF] Hess | Pre-and post-transitional behavior of the flow alignment and flow-induced phase transition in liquid crystals[END_REF][START_REF] Doi | The Theory of Polymer Dynamics[END_REF], leading to applications among which number of works should be mentionned [START_REF] Forest | Full-tensor alignment criteria for sheared nematic polymers[END_REF][START_REF] Forest | Anisotropy and dynamic ranges in effective properties of of nematic polymer nano-composites[END_REF][START_REF] Rienäcker | Chaotic and regular shear-induced orientational dynamics of nematic liquid crystals[END_REF][START_REF] Hess | Regular and chaotic orientational and rheological behaviour of liquid crystals[END_REF][START_REF] Hess | From shear-thickenning and periodic flow behavior to rheochaos in nonlinear Maxwell model fluids[END_REF].

In addition, many situations exist where the observed unsteady properties result in propagating waves or instabilities, leading to periodic motions (stick-slip like, alternating bands…) or irregular chaotic behaviour [START_REF] Herle | Alternating Vorticity Bands in a Solution of Wormlike Micelles[END_REF][START_REF] Rienäcker | Chaotic and regular shear-induced orientational dynamics of nematic liquid crystals[END_REF][START_REF] Hess | Regular and chaotic orientational and rheological behaviour of liquid crystals[END_REF][START_REF] Hess | From shear-thickenning and periodic flow behavior to rheochaos in nonlinear Maxwell model fluids[END_REF][START_REF] López-González | Shear Banding Fluctuations and Nematic Order in Wormlike Micelles[END_REF][START_REF] Cates | Rheological chaos in a scalar shear-thickening model[END_REF][START_REF] Manneville | A spatio-temporal study of rheo-oscillations in a sheared lamellar phase using ultrasound[END_REF][START_REF] Fielding | Spatiotemporal Oscillations and Rheochaos in a Simple Model of Shear Banding[END_REF].

On the other hand, it was early suggested from instabilities and hysteresis cycles observed in polymer melt extrusion that assuming a multivalued flow curve could be sufficient to generate instability [START_REF] Quemada | Unstable Flows of Concentrated Suspensions[END_REF][START_REF] Mcleish | A molecular approach to the spurt effect in polymer melt flow[END_REF]. More recently, extensive work on stability analysis has confirmed that the region where

γ σ & d d
<0 is indeed responsible of band formation [START_REF] Bautista | Stability Analysis of Shear Banding Flow with the BMP model[END_REF][START_REF] Berret | Metastable versus unstable transients at the onset of a shear-induced phase transition[END_REF][START_REF] Porte | Inhomogeneous Flows of Complex Fluids: Mechanical Instability Versus Non-Equilibrium Phase Transition[END_REF][START_REF] Dhont | A Constitutive Relation Describing the Shear Banding Transitions[END_REF][START_REF] Spenley | Nonmonotonic Constitutive Laws and the Formation of Shear-Banded Flows[END_REF].

However, most of these theoretical models fail in predicting the robustness of the stress value at the plateau, σ*, if they do not incorporate in the constitutive relation some non-local effect as diffusion terms concerning either concentration [START_REF] Olmsted | Phase Separation of Rigid-Rod Suspensions in Shear Flow[END_REF][START_REF] Fielding | Flow phase diagrams for concentration-coupled shear banding[END_REF], shear stress [42,43 ,45-47], shear rate [START_REF] Dhont | A Constitutive Relation Describing the Shear Banding Transitions[END_REF] or both shear stress and shear rate [START_REF] Wilson | Linear instability of planar shear banded flow both diffusive and non-diffusive Johnson-Segalman fluids[END_REF]. This allows to change each (unphysical) discontinuity in γ & into a "material interface", the thickness of which is of the order of the diffusion length, with γ & varying continuously within the interface from low to high values in outside bands, 1 γ& and 2 γ& . Note that considering only a coupling through concentrationdependent stress σ(φ)or viscosity η(φ)where φ is the particle volume fraction also allowed to obtain a selection of σ * as in reaction-diffusion models [START_REF] Goveas | A Minimal Model for Vorticity and Gradient Banding in Complex Fluids[END_REF].

As a contribution to phenomenological modelling of shear angular velocity under stress-banding conditions, in order to predict more precisely actual rheometrical data in a cylindrical Couette cell. This is successfully achieved in Sec.4 using the SM-model, allowing to discuss in Sec.5 an example of application to experimental data of WMS. Finally, some concluding remarks and comments on future development of this work are given in Sec.6.

-The SM-model 2.1 -Foundations of the viscosity equation

The SM-model has been deduced from minimization of viscous energy dissipation in concentrated suspensions [START_REF] Quemada | Rheology of concentrated disperse systems and minimum energy dissipation[END_REF] that led to a (Krieger-like) viscosity-volume fraction relation,

η(φ) ∼ (1-φ/φ m ) -2 .
This relation traduces a transition-like behaviour close to maximum packing of particles, φ m , with a "critical exponent" 2. For dispersions of aggregating particles, forming "Structural Units" (SU), the SU-concept led to generalize this relation to the form

η(φ eff ) ∼ (1-φ eff /φ m ) -2
, where φ eff =(1+CS)φ is the effective volume fraction of the disperse phase, S being the fraction of individual particles within all of SU, and C a factor depending on the mean compactness of SU (see [START_REF] Quemada | Rheological Modelling of Complex Fluids: I. The Concept of Effective Volume Fraction Revisited[END_REF] for further details on these definitions).

As a structural variable, S takes into account the volume of suspending fluid immobilized within SU [START_REF] Quemada | Rheological Modelling of Complex Fluids: I. The Concept of Effective Volume Fraction Revisited[END_REF][START_REF] Windhab | Fluid Immobilization -A Structure Related Key Mechanism for the Viscous Flow Behaviour of Concentrated Suspension Systems[END_REF]. Finally, shear-induced structural changes within the system are incorporated into the model assuming that S obeys a kinetic equation [START_REF]) kinetic constants κ F and κ D for formation and destruction of SU, S 0 and S ∞ being the very low-and very high-shear limits of S, respectively. The steady state solution hence depends on the ratio Γ =κ D /κ F . Again for the simplest model[END_REF]. It is assumed that, when a given shear rate is applied, a sort of order-disorder equilibrium is established, in which the forward (Brownian motion) and backward (shear-induced ordering) processes balance.

This leads to the following form for the structural viscosity

[ ] 2 S ) 1 ( 1 - ∞ - - = χ η η (1) In Eq.(1) η is a function of S=S(φ, Γ, t), where Γ = C γ γ & & is a reduced shear rate, C
γ& is a characteristic shear rate (≈ reciprocal of a structural relaxation time t C ) and

( ) 2
1 0 η η χ ∞ = ,
η ∞ and η 0 being the limiting viscosities at Γ >> 1 and Γ << 1, respectively. Note that these limits correspond to complete rupture and complete building of the structure, taking the corresponding values S ∞ =0 and S 0 =1 for simplicity (see [START_REF]) kinetic constants κ F and κ D for formation and destruction of SU, S 0 and S ∞ being the very low-and very high-shear limits of S, respectively. The steady state solution hence depends on the ratio Γ =κ D /κ F . Again for the simplest model[END_REF]).

Applying such a phenomenological model to systems like WMS is tentatively made in this work. As this model has been successfully used to describe a very large variety of concentrated dispersions [54], it is thought that some details of structural processes at the microscale would have a weak influence on rheological properties at the mesoscale. For instance, in the special case of WMS, details on reversible scission and recombination of chains and/or reptation process in the entangle state will be forgotten, keeping only the following picture : (i) the limiting structure at low shear rate is formed by large entangled groups of disordered micelles which immobilize a large amount of suspending fluid, thus having a high viscosity η 0 ; (ii) as the shear rate increases, these groups are progressively broken down with releasing some part of the entrapped fluid; (iii) finally, at high shear, both micelle alignment in the flow direction and complete release of entrapped fluid lead to a small value for η ∞ . It is expected that, despite this approach remains approximate, it benefits from a very limited number of parameters, moreover having a clear physical significance.

-Prediction of shear-banding

Under steady conditions, Γ = Cte, the shear stress variations are given by

Γ Γ χ Γ γ η γ η γ σ 2 C 1 ) (         + + = = ∞ & & & (2)
This model exhibits a non-monotonic flow curve σ(Γ) within a range of χ-values, hence leads to a multivalued shear rate for a given stress. Here we will use the variable

α = ∞ η η 0 = χ -2
for convenience. . Nonmonotonicity of flow curves Σ(Γ) occurs for α > 81 , the higher the α-value, the more pronounced the non-monotonicity.

It is interesting to observe however that the viscosity 2) is a purely monotonic function, even within the region of multivalued shear rate for a given stress.

γ γ σ γ η & & & ) ( ) ( = related to Eq.(
In terms of the reduced viscosity Η(Γ) , Fig. 3B shows behaviours that are always monotonic:

shear thinning if α > 1 (thus specially in the multivalued region) and a shear thickening if α < 1 .

-Prediction of vorticity-banding

An interesting feature of this model is the possibility to use it for modelling vorticity banding. Indeed, if the steady viscosity explicitly depends on the stress, the flow curve is defined by

) (σ η σ γ = & .
In this case, one should use a reduced shear stress Σ = σ /σ C instead of the reduced shear rate Γ in Eq.( 2), with a characteristic shear stress σ C and keeping

Γ = C γ γ & &
where

C γ& = σ C /η ∞ (i.e. still with σ C = η ∞ C γ& ), that leads to the flow curve in reduced variables, ( ) ( ) [ ] Γ Σ χ Σ Σ 2 1 + + = (3) Figure 4A displays the reduced flow curves Σ= Σ(Γ) for different values of β = 0 η η ∞ = χ 2 .
A domain of multivalued stress at given shear rate is observed for β > 81 . Corresponding curves for the reduced steady viscosity as a function of

Σ ( ) ( ) [ ] 2 1 ) ( Σ χ Σ η η Σ Η + + = = ∞ (4)
are shown on Fig. 4B. By contrast with the shear-banding case, this multivalued domain corresponds to a shear thickening behaviour. Note that in terms of shear rate the reduced viscosity Η(Γ) exhibits on Fig. 4C a multivalued domain for β > 81 as expected since a multivalued stress is then associated to a given Γ . Similar findings would be obtained for Η(Σ) in the shear banding case.

-Comparison with several minimal models

SM-model can be compared to other phenomenological ones recently proposed in the literature. In particular, the requirement for a multivalued flow curve has led to ad-hoc models, often called "minimal models", the quality of them being their simple analytical form that allows to perform stability analysis of band formation. For instance, still keeping the notation Γ for the reduced shear rate, the simplest (trinomial) dependence under steady conditions σ ~ γ & η 0 /(1+Γ 2 ) was taken for either the total shear stress (Eq.( 9) in [START_REF] Rehage | Rheological Properties of Viscoelastic Surfactant Systems[END_REF]) or its viscoelastic part only (Eqs.( 6) and ( 7) in [START_REF] Spenley | Nonmonotonic Constitutive Laws and the Formation of Shear-Banded Flows[END_REF]). A little less simple expression (also depending

on η ∞ ), σ ~ 2 2 0 1 Γ Γ η η γ + + ∞ & , (5)
has been used (Eq.( 8) in [START_REF] Rehage | Rheological Properties of Viscoelastic Surfactant Systems[END_REF] and Eq.( 33) in [START_REF] Dhont | A Constitutive Relation Describing the Shear Banding Transitions[END_REF]). The nonlinearity in Eq.( 5), resulting from the shear rate dependence in Γ 2 , leads automatically to a trinomial equation under controlled stress condition, hence to the possibility of having 3 roots. It is noteworthy that there is not other justification for selecting this kind of expression among the simplest ones. By contrast, the nonlinearity in Eq.( 2) results from the nonlinear φ-dependence of η , η ∼ (1-φ/φ m ) -2 , even (as it is considered in this paper) if the shear rate dependence of the ratio of kinetic constants remains linear (see [START_REF]) kinetic constants κ F and κ D for formation and destruction of SU, S 0 and S ∞ being the very low-and very high-shear limits of S, respectively. The steady state solution hence depends on the ratio Γ =κ D /κ F . Again for the simplest model[END_REF]). Therefore, using Eq.( 2) appears less arbitrary than using Eq.( 5).

Very similar comments can be given in the case of vorticity banding, with a constitutive equation containing an ad-hoc third order polynomial R(σ), leading to a trinomial equation in σ under controlled shear rate condition (see for instance Eqs.( 1) and ( 2) in [START_REF] Aradian | A minimal model for chaotic shear banding in shear-thickening fluids[END_REF], discarding time and non-local effects, i.e. in steady homogeneous conditions). Such an equation corresponds to the one deduced from the relation

γ σ η σ & = ) (
with η(σ) /η ∞ given by Eq.( 4).

-SM-Model characteristics under steady shear-banding conditions.

In order to predict the stress plateau associated to shear banding, the top-jumping concept have been at first introduced. Although this was done on the basis of satisfactory agreement between model predictions [START_REF] Mcleish | A molecular approach to the spurt effect in polymer melt flow[END_REF][START_REF] Spenley | Nonlinear Rheology of Wormlike Micelles[END_REF] and WMS data [START_REF] Rehage | Viscoelastic Surfactant Solutions : Model Systems for Rheological Research[END_REF], stress startup behaviour under constant γ& shows clearly that flow curve determination requires to perform measurements exceedingly slowly. Transient stress responses obtained using startup experiments showed indeed two different behaviours. Below the plateau regime ( γ& < 1 γ& ) the stress grows rapidly (in times ≤ Maxwell time τ M ) up to a stationary value and remains constant. Above 1 γ& , in the plateau region, the transient stress exhibits the same rapid increase at short times followed by a slow relaxation on very longer times τ R >> τ M , fairly described by a stretched exponential [START_REF] Berret | Rheology of Wormlike Micelles: Equilibrium Properties and Shear Banding Transition[END_REF][START_REF] Berret | Metastable versus unstable transients at the onset of a shear-induced phase transition[END_REF]. Such findings have demonstrated that reaching steady conditions requires very long delay in the plateau regime, likely due to the development of structures associated with band formation. Therefore, it appears that number of "steady" measurements, which under (up and down) imposed stress yield flow curves similar to hysteresis cycles [START_REF] Bonn | Bistability in non-Newtonian flow: Rheology of lyotropic liquid crystals[END_REF][START_REF] Ramos | Shear Melting in Lyotropic Hexagonal Phases[END_REF],

actually correspond to unsteady behaviour. Indeed, a flow-curve resembling top-jumping can be recorded by fast stress-controlled scans (roughly one data point per time ≈ 10τ M ) [START_REF] Grand | Slow Transients and Metastability in Wormlike Micelle Rheology[END_REF]. In some extent, the SM-model may predict such a unsteady behaviour, as follows.

Considering a fluid governed by Eq.( 2) with χ < 9 , one may discuss the predicted shear rate if the stress is increased fast enough (with a scanning time τ S such as τ M < τ S <<τ R ) to assume that structural changes are delayed, i.e. keeping approximatively the viscosity value of the previous step. Such an assumption is acceptable if τ S is of same order than

- = n n n γ η σ γ & &
, where the subindex n stands for the n-th data point. (i) Σ reaches progressively the apparent plateau, with a "Σ jump -value" clearly above Σ max (ii) the successive data points in response to identical stress increments dΣ correspond to a series Γ n which exhibits successive jumps with an increasing increment dΓ n = (Γ n -Γ n-1 ).

These predictions are in fair agreement with several experiments (see Figs.3a-7 in [START_REF] Grand | Slow Transients and Metastability in Wormlike Micelle Rheology[END_REF]) which exhibit both the smooth transition to an apparent plateau (well above the "true" plateau for shear banding) and data points with increasing dΓ. This confirms that stress controlled flow curves displayed on Fig. 5 are definitely nonequilibrium curves and concern successive structural states which could be considered as points of some "metastable extension" of the low shear monophasic branch.

Lastly, Fig. 6 shows that decreasing Σ from the high shear branch of the flow curve leads to a bottom-jumping-like behaviour hence, with the previous Σ-increase, to the formation of an hysteresis cycle very similar to the ones observed experimentally [START_REF] Bonn | Bistability in non-Newtonian flow: Rheology of lyotropic liquid crystals[END_REF][START_REF] Ramos | Shear Melting in Lyotropic Hexagonal Phases[END_REF].

As yet discussed in Sec.1, top or bottom jumpings thus do not give the good criterion for selecting σ * . This can be only achieved by using theories which involve constitutive equations with nonlocal terms. Therefore, we do not proceed to further development of the SM-model before introduction of such non-local effects. In the next sections, the discussion will be limited to the applicability of this model to the prediction of some data.

4-Equations of Couette rheometry

After a brief review of the problem related to the determination of the shear rate function in monophasic flows, a special attention is devoted to the calculation of the torque vs. angular velocity under stress-banding conditions, in order to predict more precisely actual rheometrical data in a cylindrical Couette cell. This is successfully achieved using the SMmodel.

-Monophasic flows

Determining the flow curve ) (γ σ & of non-Newtonian fluids requires the knowledge of both the shear stress σ and the shear rate γ & at one place in the flow domain of the cell. In rheometry between concentric cylinders (discarding end effects), it is worth to observe that velocity vector v and shear stress tensor σ in cylindrical coordinates [r, θ, z] reduce [61] to their components v θ = v(r,t) and σ rθ = σ(r,t), leading to the balance condition σ(r,t)

= ) t ,
∫ = R R dr r r Ω κ γ ) ( & , (7)
where R and κR are the radii of outer and inner cylinder, respectively (0 < κ < 1), and L is the cylinders height. While σ(r) is obtained straightforwardly from Eq. ( 1), extracting

) (r γ &
requires inverting the integral of Eq.( 7), the solution of which is not unique due to the scattering present in experimental data. This inverse problem is also designated ill-posed in the literature [START_REF] Friedrich | New Ill-Posed Problems in Rheology[END_REF][START_REF] Berli | A Procedure to Determine the Viscosity Function from Experimental Data of Capillary Flow[END_REF]. Indeed, direct estimations of γ & can be achieved only when the gap between inner and outer cylinders is very small, say 0.99 < κ < 1. Under these circumstances,

) (r γ &
is considered nearly uniform throughout the flow domain, and thus Eq.( 7) can be replaced by ) 1

( κ κ Ω γ - ≈ ng &
, which is known as the narrow-gap solution for Couette viscometry [62,54].

For wider gap widths, an approach commonly used in practice consists in estimating the shear rate at the inner wall,

) ( R R κ γ γ κ & & =
, by using the following asymptotic expression [START_REF] Krieger | Direct Determination of the Flow Curves of Non-Newtonian Fluids. II. Shearing Rate in the Concentric Cylinder Viscometer[END_REF],

... ln d d 3 ) / 1 ln( ln d d ) / 1 ln( 2 R 2 R R + + + = κ κ κ σ Ω κ σ Ω κ Ω γ& , (8) where ) 2 ( 2 2
L R M R πκ σ κ =
is the shear stress at the inner cylinder. This equation is recommended in rheology books (for instance, [54]) to be used for 0.5 < κ < 1. It is relevant to mention that the error in terminating the series at the second term is lower than [ ] 2

) / 1 ln( κ .

As an interesting alternative to estimate R κ γ & , finite difference formulas that avoid the numerical differentiation of raw data were recently proposed [START_REF] De Hoog | Simple and Accurate Formulas for Flow-Curve Recovery from Couette Rheometer Data[END_REF]. Although the observed robust selection of the stress plateau σ = σ * is only obtained with theories including non-local (diffusive) terms in the constitutive model, we intend to show that it is possible to limit the modelling of stress controlled measurements to these outside bands, considering the plateau value σ* as a model variable to be adjusted together with the SM-model parameters, η ∞ , η 0 and t C ≡ C γ& -1 .

-Shear-banding flows

In Couette cells, the shear stress is maximum at the inner cylinder ( R κ σ ), and decreases smoothly as

2 ) / ( ) ( r R r R κ σ σ κ =
, to give the minimum value at the outer cylinder

( R R κ σ κ σ 2 =
). The simple (generally accepted) scenario considered here is that, when the In principle, these assumptions agree with hypothesis and experimental data reported in the literature ( [START_REF] Cappelaere | Rheology, Birefringence, and Small-Angle Neutron Scattering in a Charged Micellar System: Evidence of a Shear-Induced Phase Transition[END_REF][START_REF] Salmon | Velocity Profiles in Shear-Banding Wormlike Micelles[END_REF][START_REF] Spenley | Nonlinear Rheology of Wormlike Micelles[END_REF][START_REF] Radulescu | Matched Asymptotic Solutions for the Steady State Banded Flow of the Diffusive Johnson-Segalman Model in Various Geometries[END_REF] for instance). However, in relation with the interface (ii), very recent theoretical [START_REF] Fielding | Nonlinear dynamics of a shear banding interface[END_REF] and experimental [START_REF] Lerouge | Interface Instability in Shear-Banding Flow[END_REF] works suggest the existence of a region of instability between the bands in the vorticity direction. As a first approximation, here we assume a sufficiently narrow and flat interface between the bands. Under these conditions, and considering the SM-model, Eq.( 2), defining ) ( γ σ & as that plotted in Fig. 2, the expressions of Ω as a function of σ κR in the different stress domains are given as follows (see also [START_REF] Radulescu | Matched Asymptotic Solutions for the Steady State Banded Flow of the Diffusive Johnson-Segalman Model in Various Geometries[END_REF]), *

σ σ κ ≤ R , ∫ < - = R R d Ω σ σ κ σ σ σ γ ) ( 2 1 & ; (9a) R R σ σ σ κ ≥ ≥ * ,       + - = ∫ ∫ > < * * ) ( ) ( 2
1 σ σ σ σ κ σ σ σ γ σ σ σ γ R R d d Ω & & ; (9a) * σ σ ≥ R , ∫ > - = R R d Ω σ σ κ σ σ σ γ ) ( 2 1 & . ( 9c

-Calculations proposed for shear-banding flows

-Direct and inverse calculations

For a given fluid, if the parameters (η 0 , η ∞ , t c ) entering the SM-model are known, Eqs.

(2) and(9a)-(9c) predict the values Ω vs. σ (or M) to be obtained in an experiment in which the requirements to accomplish a viscometric flow are satisfied (steady state, no-slip at the walls, end effects negligible, isothermal flow), as well as the assumptions (i-iii) made above. This possibility, designated 'direct calculation', will be illustrated in detail in Section 5.

On the other hand, a more challenging problem is the inverse calculation: determining the values of the parameters (η 0 , η ∞ , t c ) from the curve of raw data Ω vs. σ (or M), and then using them to plot the flow curve ) (γ σ & in the appropriate range of shear rates. The implementation of this task requires further efforts, because the minimization problem involved cannot be solved with standard mathematical software (for instance, [START_REF] Berli | A Procedure to Determine the Viscosity Function from Experimental Data of Capillary Flow[END_REF]). for a given set of known parameters (η 0 , η ∞ , t c ). This task is carried out through a Newton-Raphson subroutine [START_REF] Carnahan | Applied Numerical Methods[END_REF], for around 10 4 discrete values of σ in the range

-Numerical procedure for the direct calculation

R R κ σ σ σ ≤ ≤
. Then Eqs.(9a)-(9c) are integrated numerically by using the trapezoidal rule, also with 10 4 discrete intervals. This is carried out for an initial, tentative value * σ (normally inferred from the experimental curve), which is then adjusted to provide the best representation of data.

One should underline that the present calculation takes into account the variations of

) (σ γ < & and ) (σ γ >

&

(i.e. spatial variations) in each band. Hence, it differs notably from current assumptions made in several works [START_REF] Cappelaere | Rheology, Birefringence, and Small-Angle Neutron Scattering in a Charged Micellar System: Evidence of a Shear-Induced Phase Transition[END_REF][START_REF] Porte | Inhomogeneous Flows of Complex Fluids: Mechanical Instability Versus Non-Equilibrium Phase Transition[END_REF][START_REF] Cates | The Rheology of Entangled Polymers at Very High Shear Rates[END_REF][START_REF] Georgiou | On the Stability of the Simple Shear Flow of a Johnson-Segalman Fluid[END_REF][START_REF] Drappier | Plateaux de contrainte et bandes de cisaillement dans les fluides complexes[END_REF]. In these studies, < being the gap width, and γ& being the "measured" shear rate. Such constant values are not observed in measured velocity profiles (see, for example, Fig. 2 in [START_REF] Cappelaere | Rheology, Birefringence, and Small-Angle Neutron Scattering in a Charged Micellar System: Evidence of a Shear-Induced Phase Transition[END_REF]). Moreover, it is evident that the larger the gap thickness e, the higher the differences from constant values of γ & .

-Examples of application

This section illustrates the applicability of the SN-model to represent rheometric data of shear-banding fluids. The first step consists in confronting Eq.( 2) to data σ vs. γ & from real systems. For this purpose, experimental data of WMS published in the literature are considered: CTAB [START_REF] Cappelaere | Rheology, Birefringence, and Small-Angle Neutron Scattering in a Charged Micellar System: Evidence of a Shear-Induced Phase Transition[END_REF] and CPCl-NaSal [START_REF] Salmon | Velocity Profiles in Shear-Banding Wormlike Micelles[END_REF]. These data, which were reported as R κ σ vs. ng γ & , are presented in Fig. 7 as

R κ σ vs. κ κ γ Ω / ng ) 1 ( - = &
, by using the corresponding values of κ (0.94 and 0.96, respectively).

Applying Eq.( 2 In order to cross-check these results, and having the parameters η 0 , η ∞ , and t c that characterize the fluid, we finally carried out the direct calculation, defined in Sec.

-Concluding remarks and perspectives

A phenomenological (structural) model has been introduced in order to improve rheological modelling of complex fluids under shear banding conditions. In contrast with the majority of other minimal models already developed in the literature for shear banding description, this model has a "structural origin" (a kinetic description which accounts for shear induced changes of the internal structure of the material) which seems to bring forwards γ& to be adjusted by data fitting.

A crucial aspect in complex fluids like WMS, and more generally systems exhibiting shear banding phenomena, is the modelling of unsteady rheometric data, specially stress relaxation after the sudden inception of a given shear rate. A stress overshoot is observed, followed by dumped oscillations and either or not Maxwellian relaxation to the steady state [START_REF] Lerouge | Correlations between Rheological and Optical Properties of a Micellar Solution under Shear Banding Flow[END_REF][START_REF] Berret | Rheology of Wormlike Micelles: Equilibrium Properties and Shear Banding Transition[END_REF][START_REF] Berret | Metastable versus unstable transients at the onset of a shear-induced phase transition[END_REF]. Theoretical descriptions of these experiments are rather demanding, as different time-dependent phenomena are involved, mainly fluid microstructure relaxation and fluid viscoelasticity, even the transient characteristics of the apparatus in some studies [START_REF] Talbi-Boucenna | Transient Measurements in Rate Controlled Rheometers[END_REF]. The structural model from which Eq.( 2) derives also includes shear induced kinetic processes involved in structural changes (see [START_REF]) kinetic constants κ F and κ D for formation and destruction of SU, S 0 and S ∞ being the very low-and very high-shear limits of S, respectively. The steady state solution hence depends on the ratio Γ =κ D /κ F . Again for the simplest model[END_REF]). A nonlinear Maxwell model based on Eq.( 2) allowed to interpret unsteady curves on the basis of such kinetic processes [54,[START_REF] Quemada | Rheological modelling of complex fluids: IV: Thixotropic and "thixoelastic" behaviour. Start-up and stress relaxation, creep tests and hysteresis cycles[END_REF], hence to predict transient responses like those observed in WMS. This work is in progress.

Future work will be devoted to investigate if the SM-model predictions after including non-local terms will make possible selecting the σ*-value and further would provide some new physical argument in the debate on the origin of shear banding (i.e. formation and robustness of the stress plateau) in favour of either a mechanical instability or a shear-induced phase transition [START_REF] Porte | Inhomogeneous Flows of Complex Fluids: Mechanical Instability Versus Non-Equilibrium Phase Transition[END_REF]. Moreover, these model predictions will be aimed to be compared with those from Johnson-Segalman model, which is nowadays the most employed constitutive relation to describe shear-banding fluids [START_REF] Fielding | Spatiotemporal Oscillations and Rheochaos in a Simple Model of Shear Banding[END_REF][START_REF] Olmsted | Phase Separation of Rigid-Rod Suspensions in Shear Flow[END_REF][START_REF] Fielding | Flow phase diagrams for concentration-coupled shear banding[END_REF][START_REF] Fielding | Linear instability of planar shear banded flow[END_REF][START_REF] Fielding | Nonlinear dynamics of a shear banding interface[END_REF][START_REF] Aradian | A minimal model for chaotic shear banding in shear-thickening fluids[END_REF][START_REF] Radulescu | Matched Asymptotic Solutions for the Steady State Banded Flow of the Diffusive Johnson-Segalman Model in Various Geometries[END_REF][START_REF] Georgiou | On the Stability of the Simple Shear Flow of a Johnson-Segalman Fluid[END_REF], as well as the so-called BMP model [START_REF] Bautista | Stability Analysis of Shear Banding Flow with the BMP model[END_REF].

Figure captions

  Fig.3A displays these variations in terms of a reduced stress Σ = σ /σ C vs. Γ (taking σ C = η ∞ C γ& for the sake of simplicity, without reducing the generality of the results). As yet underlined above, Eq.(2) defines a (local) reduced steady viscosity Γ, α) as an explicit function (i) of shear rate γ& through the reduced variable Γ , (ii) of volume fraction through the φ-dependent ratio of limiting viscosities, α = ∞ η η 0

Fig.5A displays

  Fig.5A displays the resulting variations of Σ vs. Γ (same variables yet introduced in the previous section) under controlled stress conditions. A top-jumping-like behaviour is clearly observed after reaching the flow curve maximum Σ max . However, as shown on Fig.5B, details of variations Σ(Γ) close to the maximum of the flow curve show two main characteristics :

  Shear-banding flows involve an additional complexity: abrupt changes in the function ) (r γ & arise when two or more phases coexist in the flow domain of the cell. As described in the introduction, different theoretical models aimed to predict the occurrence of shearbanding, with a general consensus that the underlying constitutive curve of shear-banding fluids has the form shown in Fig.2.

.

  value * σ , a new phase with lower flow resistance develops from the inner cylinder, and consequently R κ γ & jumps from 1 γ & to 2 γ & (Fig. 2). Further increase in R κ σ results in the growth of this new phase, with a thickness enhancement of the associated band, up to complete development of the same phase in the whole gap, as R For the purpose of describing our rheometric problem, we start with the following assumptions (i) shear-bands are stable and can coexist in steady state conditions;(ii) the bands present a localized, flat and sufficiently narrow interface at a certain * u(r) are continuous at the interface.

  respectively. To handle the inverse problem given by Eqs.(9a)-(9c), we need to introduce the constitutive equation ) (γ σ & characterizing the fluid under study.

γ & and > γ & are taken as the limits 1 γ & and 2 γ

 12 & of the plateau σ* , hence are the constants entering the "lever rule" e 1 and e 2 are the respective band thickness, R R e κ -=

  ) firstly requires the estimation of accurate values R κ σ vs. R κ γ& from data R κ σ vs. Ω . Thus we selected intervals of experimental data of Fig.7 where the fluid is monophasic only. The low and high shear zones were analysed independently one another to find R κ γ & by means of Eq.(6). Results are presented in Fig.8 (symbols). Afterwards, Eq.(2) is adjusted to data R κ σ vs. R κ γ& (full lines in Fig.8), leading to the following parameter values: CPCl/NaSal : η 0 = 63.6 Pa.s , η ∞ = 0.64 Pa.s , t c = 32 ms ; CTAB : η 0 = 7.8 Pa.s , η ∞ = 25 mPa.s , t c = 2.7 ms. It is observed that the model describes satisfactorily the flow curve in the full range of shear rates, by predicting an intermediate multivalued zone.

4 . 3 .

 43 That is, Eqs.(9a-c)-(2) were calculated numerically, and the resulting function ) ( experimental curve. This is actually done in Fig.9, where full lines are the numerical predictions with the adjusted values of * σ indicated in the figure caption. A remarkable agreement is observed in the full range of experimental data. Furthermore, it is worthy of note that the values of * σ obtained are in close agreement with those previously reported[START_REF] Cappelaere | Rheology, Birefringence, and Small-Angle Neutron Scattering in a Charged Micellar System: Evidence of a Shear-Induced Phase Transition[END_REF][START_REF] Salmon | Velocity Profiles in Shear-Banding Wormlike Micelles[END_REF].

a

  better physical basis. Simplest assumptions on structure kinetics lead to a Strucural Minimal (SM-) model. Existence of a flow curve, multivalued in shear rate -a necessary requirement for the presence of shear banding-, has been found depending only on the ratio of the limiting viscosities at very high and very low shear. Under controlled (increasing) stress conditions, the SM-model predicts a top-jumping-like behaviour with a smooth transition to an apparent plateau. The later cannot represent the true plateau associated to shear banding since such predictions concern nonequilibrium responses as those that could result from too fast measurements. Similarly, for decreasing stress, a bottom-jumping is predicted, forming with top-jumping a loop which resembles a hysteresis cycle. However, such findings have not been further discussed in the paper since many works have demonstrated a lack of stability in the vicinity of the stress plateau if non-local effects are absent from the model used, as it is the case in this work. Prediction of given rheometric data of shear-banding fluids has been a test for this constitutive model of the flow curve. It may be remarked that the success of calculations suggested relies on both the introduction of a suitable model for the fluid, and the adequate computation of the shear rate function in the rheometric cell. Indeed, determining the true shear rate attained in Couette flows is a non-trivial task, except if the gap between cylinders is very small. When fluids presenting strong shear-thinning behaviour are studied, an inappropriate estimation of the shear rate leads to considerable errors, notably when two phases coexist in the flow domain. In this work, although the absence of non-local terms in the SM-model impedes the stress plateau value σ* to be determined, successfully testing of the model has been nevertheless obtained by adding σ* to the three model variables η ∞ , η 0 and C

Figure 1 :

 1 Figure 1: Curve of rheometric data (arbitrary drawing) typically found in WMS studied by Couette flow. The insets are highly schematic representations of the solution structure in different shear rates zones. At a given shear stress, a new fluid phase (band) develops as micelles align in the flow direction.

Figure 2 :

 2 Figure 2:Shear stress a function of shear rate (arbitrary units) for shear-banding fluids.

Figure 3 :

 3 Figure 3:A -Reduced stress Σ vs. reduced shear-rate Γ from Eq.(2) for different values of the limiting viscosity ratio α =

Figure 4 :

 4 Figure 4: A -Reduced stress Σ vs. reduced shear-rate Γ , defined in Eq.(3), for different values of the limiting viscosities ratio β =

Figure 5 :

 5 Figure 5: A-Reduced stress vs. reduced shear rate under controlled stress conditions, showing a topjumping -like behaviour (△) (Stress increment : dΣ = 4.10 -2 ), in comparaison with the

Figure 6 :

 6 Figure 6:Idem Fig.5, showing top-and bottom-jumpings, leading to a hysteresis cycle.

Figure 7 :

 7 Figure 7: Angular velocity as a function of controlled shear stress for different WMS. Symbols represent two examples of experimental data reported in the literature: CTAB [6], and CPCl/NaSal [7].

Figure 8 :

 8 Figure 8:Shear stress as a function of shear rate, corresponding to the systems presented in Fig.7.Symbols are the values obtained from experimental data ) (Ω σ κR

Figure 9 :

 9 Figure 9: Shear stress as a function of angular velocity for different WMS. Symbols represent two examples of experimental data reported in the literature: CTAB [6], and CPCl/NaSal [7]. Full lines are the predictions of Eqs.(9a-c)-(2), with the values of η 0 , η ∞ , and t c reported in the text. In addition, * σ = 62 Pa for CPCl/NaSal, and * σ = 32.5 Pa for CTAB.

  Figure 1 Quemada and Berli 2008

  Figure 6 Quemada and Berli 2008

& & &

0,14 0,16 0,18 0,2 0,22

0,00001 0,001 0,1 10

Acknowledgements: C. B. acknowledges CONICET and SEPCYT-FONCyT, Argentina, for the financial aid received.

imposed stress (up) imposed stress (down) imposed shear rate


Esteban Carvalho

Pierre Susbielle

Ahmad Hably

Jilles S Dibangoye

Nicolas Marchand

Estéban Carvalho

Neural Enhanced Control for Quadrotor Linear Behavior Fitting

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

I. INTRODUCTION

In the last decade, unmanned aerial vehicles (UAVs) have become predominant worldwide. From leisure film making to public transport, surveillance, monitoring or even geographic mapping, their area of application grows continuously. The flight navigation is usually done manually, which requires a minimum of training to pilot UAVs. A key component to make UAVs available for all is an automatic flight control reactive to unexpected situations and environments.

A functional autopilot consists of two elements: a trajectory generation [START_REF] Richter | Polynomial trajectory planning for aggressive quadrotor flight in dense indoor environments[END_REF][START_REF] Brescianini | Computationally efficient trajectory generation for fully actuated multirotor vehicles[END_REF] and a controller [START_REF] Bangura | Nonlinear dynamic modeling for high performance control of a quadrotor[END_REF][START_REF] Faessler | Differential flatness of quadrotor dynamics subject to rotor drag for accurate tracking of high-speed trajectories[END_REF]. On the first hand, a fast trajectory generation algorithm is crucial to maintain good flight performances. It should be able to face unpredicted obstacles and respond in real time. As a result, [START_REF] Richter | Polynomial trajectory planning for aggressive quadrotor flight in dense indoor environments[END_REF][START_REF] Brescianini | Computationally efficient trajectory generation for fully actuated multirotor vehicles[END_REF] propose a polynomial path generator that aims rapid and feasible trajectory generation, with lowcost computation and complexity. On the other hand, the controller tracks the desired trajectory. Research of [START_REF] Bangura | Nonlinear dynamic modeling for high performance control of a quadrotor[END_REF] introduce an accurate quadrotor model that fits the system dynamic and thus improves flight performances. However, it induces at the same time a nonlinear controller which leads complex to tune algorithm as well as computation consuming systems. Although a linear model is easier to handle, a popular approach consists in linearizing a nonlinear model to avoid complexity [START_REF] Martins | Linear quadratic regulator for trajectory tracking of a quadrotor[END_REF].

*This work has been partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01) and ROBOTEX 2.0 (Grants ROBOTEX ANR-10-EQPX-44-01 and TIRREX ANR-21-ESRE-0015) funded by the French program Investissements d'Avenir.

1 Univ. Grenoble Alpes, CNRS, Grenoble INP, GIPSA-Lab, F-38000, Grenoble, France; email: firstname.lastname@gipsa-lab.fr 2 Univ. Lyon, INSA Lyon, CITI, F-69621, Villeurbanne, France; email: firstname.lastname@inria.fr With the advances in artificial intelligence, modelling and control are currently evolving. As presented in survey [START_REF] Hunt | Neural networks for control systems -A survey[END_REF], combining machine learning techniques with traditional control approaches is improving systems' performances. Data-driven methods achieve suitable modelling of errors, uncertainties and perturbations. For instance, [START_REF] Drews | Aggressive deep driving: Combining convolutional neural networks and model predictive control[END_REF] highlights the benefits of Convolutional Neural Networks (CNN) in a fast driving car control. Works [START_REF] Li | Deep neural networks for improved, impromptu trajectory tracking of quadrotors[END_REF]- [START_REF] Torrente | Data-driven MPC for quadrotors[END_REF] focus on UAVs hybrid approaches. In [START_REF] Torrente | Data-driven MPC for quadrotors[END_REF], Gaussian processes are introduced to model disturbances which are then combined to a predictive model control. This greatly reduces tracking errors at high speeds. Authors of [START_REF] Li | Deep neural networks for improved, impromptu trajectory tracking of quadrotors[END_REF] propose a Deep Neural Network (DNN) to generate a reference trajectory. It shows a gain in tracking with smaller errors while keeping classical PID controller. Finally, a predictive control for racing drone is presented in [START_REF] Kaufmann | Beauty and the beast: Optimal methods meet learning for drone racing[END_REF]. The vehicle flies through moving gates. Gates' position is predicted by a CNN and sent to the control. With this algorithm, the team managed to beat the second team in half the time.

One may use neural networks to model system dynamics, as shown in [START_REF] Mohajerin | Modelling a quadrotor vehicle using a modular deep recurrent neural network[END_REF], and then apply the learned model to adapt a control law. [START_REF] Abbeel | Autonomous helicopter aerobatics through apprenticeship learning[END_REF][START_REF] Punjani | Deep learning helicopter dynamics models[END_REF] have shown the relevance of using learning on modelling complex aerodynamics phenomena usually neglected, such as vibrations, for helicopter flight control. Based on those advances, [START_REF] Bansal | Learning quadrotor dynamics using neural network for flight control[END_REF] came up with a shallow neural network (NN) that model UAVs dynamics. A more recent study [START_REF] Shi | Neural lander: Stable drone landing control using learned dynamics[END_REF] applies a DNN to model ground effects and proposes a DNN-based controller that enhances the landing behavior, achieving impressive results.

The problem of getting an efficient controller that fits some given expected behavior can also be addressed by using reinforcement learning (RL) methods. RL algorithms for drone control is developing a lot, see [START_REF] Hwangbo | Control of a quadrotor with reinforcement learning[END_REF][START_REF] Lambert | Low-level control of a quadrotor with deep model-based reinforcement learning[END_REF]. However, these contributions fall short in that they need a huge amount of data. In addition, the data collected must be relevant to learn a controller and thus requires a large exploration process, usually done randomly. Hence for high order systems like UAVs, it requires a lot of data and time. Secondly, those unsupervised method do not provide stability guaranty for the system which is very problematic. To overcome the data collection and time problem, as well as the stability issue, we propose an hybrid control approach based on a linear control and a neural network that describes nonlinear error dynamics. Our approach allows to easily setup a controller that will fit expected behavior. In order to reduce the design complexity and to preserve stability, a linear cascaded control is used. A Deep Neural Network is added for fast performance enhancement using collected data. It will learn the error between the linear model and the real system dynamics, both controlled with the linear law. Fitting an a priori given linear behaviour is especially of interest for obstacle avoidance. Indeed, it is far easier to predict the error between a reference trajectory and the real one when the system is linear. It also eases a fast trajectory generation. The DNN will not only model a complex and a nonlinear dynamic error, but also includes parametric errors. The error model identified by the DNN is then injected in the real time command so that the quadrotor behavior fit to the linear expected dynamics. Proposed methodology is summarized in fig. 1 Main contributions. In contrast to prior research, the contributions of this paper are threefold:

• Learning the error between the real system dynamics and a linear dynamical model. • We consider a linear control applied to a complex system, avoiding nonlinear controls, and correcting the resulting error to refit the linear expected quadrotor dynamics. • We show good results in simulation and in real time experimental flights and propose an open source ROS package. The paper is organized as follows. Section II defines the quadrotor dynamics and its linearized model used for control design. Section III presents a neural architecture to model the error dynamics. Proposed linear cascaded controller and DNN correction law are given in Section IV. We expose simulation and experimental results in Section V, including neural networks training and controller's tests. Finally, conclusion and future work perspectives are given in Section VI.

II. QUADROTOR DYNAMICS

Quadrotors are aerial vehicles built with two counter-rotating rotors located on a rigid cross-shaped body. Translations on z axis are achieved by decreasing or accelerating the four rotor speeds at constant thrust for each of them. Translations on the x or y axis are obtained by adjusting the pitch or roll by applying more thrust to the two adjacent rotors and less to the two opposite ones. Finally, yaw is achieved by applying more thrust to the rotors rotating in the desired direction and less to the other two. In this section, the quadrotor's Newton-Euler dynamic and the linear model used for control purposes in Section IV are presented.

A. Quadrotor Dynamics

As a starting point, {I} defines the inertial frame described by its units vectors { i 1 , i 2 , i 3 } and {B} represents the body fixed frame described by its units vectors { b 1 , b 2 , b 3 }. R ∈ SO(3) defines the orientation of the body frame {B} in {I} using Z-Y-X Euler formalism(see Fig. 2): where c • := cos(•) and s • := sin(•).

R = R z (ψ)R y (θ)R x (φ) = c θ c ψ s φ s θ c ψ -c φ s ψ c φ s θ c ψ +s φ s ψ c θ s ψ s φ s θ s ψ +c φ c ψ c φ s θ s ψ -s φ c ψ -s θ s φ c θ c φ c θ (1)
Let ξ := (x, y, z) T ∈ {I} be the position of the drone center of mass, v := (v x , v y , v z ) T ∈ {I} its linear velocity expressed in {I}, Ω := (p, q, r) T the angular velocity expressed in {B} and ζ := (φ, θ, ψ) denotes Euler angles. The total mass of the vehicle is given by m and the inertia matrix J assumed diagonal due to the symmetry of the quadrotor. The quadrotor rigid body dynamics are given by Newton-Euler equations [START_REF] Mahony | Multirotor aerial vehicles: Modeling, estimation, and control of quadrotor[END_REF]:

           ξ = v, m v = -mg i 3 + RT b 3 + δ v nd , Ṙ = RΩ × , J Ω = -Ω × J Ω + Γ + δ Ω nd , (2a) (2b)
(2c) ( 2d
)
where T is the total thrust generated by the four rotors, Ω × is a skew-symmetric matrix associated to Ω that verifies for all matrix M ∈ M 3 (R), M × Ω = M Ω × , and Γ represents the torques applied to the rotors, and δ v nd and δ Ω nd the unmodelled and unknown dynamics on linear and angular acceleration respectively. It will be referred as error dynamics in the following. These error dynamics include for instance the gyroscopic effects, ground effects, blade flapping or battery discharge.

In the following sections, we assume that the angular loop of the quadrotor is controlled, with a sufficient speed for the other cascaded control loops. As a result, the model to be controlled is given by:

       ξ = v, m v = -mg i 3 + RT b 3 + δ v nd , ζ = W -1 (φ, θ, ψ)Ω, (3a) (3b) (3c)
where:

W -1 (φ, θ, ψ) = 1 s φ t θ c φ t θ 0 c φ -s φ 0 s φ /c θ c φ /c θ . ( 4
)
Summarizing previous equations, the quadrotor dynamics can be expressed as ( 5):

Ẋ = f X4 (X, U ) + δ nd , ( 5
)
where X is the state vector defined by

X := ξ v ζ T ,
the control vector U := T p q r T , the error

dynamics vector δ nd := 0 T 3 δ v nd m T 0 T 3 T
, with 0 3 := 0 0 0 T and:

f X4 (X, U ) =   v -g i 3 + RT b3 m T r (φ, θ, ψ)Ω   . (6)

B. Quadrotor Linear Dynamics Model

The proposed control law is derived from a linearized model around given equilibrium points. The quadrotor linear model usually used for control design is obtained from hover conditions. Equilibrium points are given by ( 7) and ( 8):

X eq = x e y e z e 0 3 0 3 T , (7)
U eq = mg 0 3 T . ( 8
)
We can then define variation variables X := X -X eq and Ũ := U -U eq . Assuming δ nd is negligible, δ v nd ≈ 0 3 , the linear model obtained using first order Taylor expansion is given by:

Ẋ = 03 I3 03 03 03 A2,3 03 03 03 X + 034 B2 B3 Ũ , Ẋ = A X + B Ũ , (9)
where:

A = ∂f X4 ∂X (X, U ) (Xeq,Ueq) , B = ∂f X4 ∂U (X, U ) (Xeq,Ueq)
, 0 3 refers to the square null matrix of order 3 and I 3 the identity matrix of order 3. 0 34 is a null matrix composed of 3 lines and 4 columns. A 2,3 , B 2 and B 3 matrices are expressed as follow:

A 2,3 = 0 g 0 -g 0 0 0 0 0 , B 2 = 0 0 0 0 0 0 0 0 m -1 0 0 0 , B 3 = 0 1 0 0 0 0 1 0 0 0 0 1 .

III. DEEP NEURAL NETWORK STRUCTURE

This section presents the proposed structure for the Deep Neural Network. Networks are used to describe the error dynamics defined by eq. ( 10), corresponding to the unknown dynamics between the real drone dynamics and the linear acceleration dynamics (δ v e ).

δ v e (ζ, T ) = v -(A 2,3 ζ + B 2 Ũ ), (10)
According to the universal approximation theorem, one can build a feed-forward network composed of only one hidden layer with a finite number of units to approximate our continuous functions δ v e given by [START_REF] Torrente | Data-driven MPC for quadrotors[END_REF]. Here, it is proposed to use a deep neural network to better fit the data. Based on the nature of the problem, here regression of temporal data, the use of this type of network is well suited.

The activation function used for all hidden layers is the rectified linear unit (ReLU) defined element-wise by h(•) = max(•, 0). Weights and bias are defined as follows:

W v,1 ∈ R #Ev,Nv , W 2 ∈ R Nu,Nu , W 3 ∈ R Nu,#Sv M 1 ∈ R Nu , M 2 ∈ R Nu and M 3 ∈ R #Sv . N u
is the number of units of the hidden layer of the network and using the following notation #A denotes the cardinal of A (fig. 3). E DN N denotes the neural network inputs. Our neural network structure can then be expressed as [START_REF] Mohajerin | Modelling a quadrotor vehicle using a modular deep recurrent neural network[END_REF]:

δv nd (E DN N ) = W 3 T h(W 2 T h(W 1 E DN N + M 1 ) + M 2 ) + M 3 . ( 11
)
We choose N l = 2 for the number of hidden layers and N u = 64 for the number of hidden layers units, as it is sufficient to correctly model the error dynamics and to perform real-time computation. A quick ablation study revealed that for one layer DNN, the network has difficulty to converge and for three layers it does not bring better performances while increasing the learning process time. To learn error accelerations dynamics one aims at finding weights and bias that minimize the mean squared error (MSE) between predictions and observations. Predictions are expected values from the neural network model and observations are the values obtained via measurements. Accelerations measurements are collected from the inertial measurement unit (IMU) of the quadrotor. Take ϑ := (W 1 , W 2 , W 3 , M 1 , M 2 , M 3 ) the parameter to be optimized. Then, training the neural network [START_REF] Mohajerin | Modelling a quadrotor vehicle using a modular deep recurrent neural network[END_REF] aims at solving the following optimization problem [START_REF] Abbeel | Autonomous helicopter aerobatics through apprenticeship learning[END_REF] over the training data-set:

min ϑ L(ϑ) = min ϑ 1 N N k=1 δv e,k (ϑ) -δ v e,k 2 , ( 12
)
where N is the number of elements in the data-set and δ v e,k are error dynamics computed by taking measurements outputs and commands inputs. Here L(ϑ) is the objective function to be minimized, referred as loss function. This problem can be typically solved using traditional gradient descent methods using back-propagated gradients optimizers.

IV. QUADROTOR CONTROL

In this section, we describe the cascaded controllers designed to track a given trajectory reference (See fig. 1, step 1). A linear cascaded architecture is used for tuning simplicity and ease of implementation. This control will be next enhanced with a deep neural network feedforward correction term. The DNN correction term aims to get closer to the linear expected behavior. This term is directly designed from the previously built neural network. Thus, one is able to rectify unmodelled dynamics and linearization errors, that are learned dynamics with training data.

A. Linear cascaded controller

The proposed control law is divided into two linear cascaded controllers: the attitude and the position/speed control loops.

Attitude controller: Using (9), it comes that: ζ = Ω. We can directly choose:

Ω = -K ζ φ θ ψ - φ θ ψ ref , (13)
     ẋ1 = x 2 , ẋ2 = 0 g 0 -g 0 0 0 0 m -1 u ζ = A ζ u ζ . (14)
Defining η 1 := x 1 -x 1ref and η 2 := η1 := ẋ1 -ẋ1ref = x 2x 2ref , ( 14) can be rewritten as:

η = 0 3 I 3 0 3 0 3 η + 0 3 A ζ (u ζ -A -1 ζ ẋ2ref ). (15)
A = 03 I3 03 03 , B = 03

Aζ

For [START_REF] Shi | Neural lander: Stable drone landing control using learned dynamics[END_REF], we propose to use again a linear quadratic regulator:

u * ζ = -K η η + A ζ ẋ2 ref , (16)
where K η is chosen such that we minimize a quadratic cost expressed as follows:

J(u ζ ) = +∞ 0 (η T (t) Q η(t) + u T ζ (t) R u ζ (t))dt.
Weighting matrices Q and R are tuned to reach desired tracking performances.

B. DNN Correction Controller

According to eq. (2b), the real system dynamics can be written as follows:

η = Aη + B(u ζ -A -1 ζ ẋ2ref + A ζ -1 δ v e ). (17)
The DNN controller is designed as follows:

u ζ = u * ζ + ∆u ζ , (18)
where u * ζ is the linear quadratic control ( 16) and ∆u ζ is the feedforward correction term based on the neural network. The closed loop form of ( 17) using ( 18) is:

η = (A -BK η )η
Close loop linear behavior

+ B(∆u ζ + A -1 ζ δ v e ).
Term to minimize [START_REF]4-autopilot software[END_REF] According to [START_REF]4-autopilot software[END_REF], the closed loop behavior is the sum of two terms. The first one describes the desired linear behavior, which is fully controlled. The second term needs to be canceled or at least minimized to get closer to the expected linear behavior. This procedure outlines a classical minimization problem [START_REF] Chollet | Keras[END_REF]:

min ∆u ζ ∆u ζ + A -1 ζ δ v e (∆u ζ ). (20)
Replacing error dynamics by their approximations from the DNN we get:

min ∆u ζ ∆u ζ + A -1 ζ δv nd (∆u ζ ). (21)
As δv nd depends on ∆u ζ according to [START_REF] Dozat | Incorporating Nesterov Momentum into Adam[END_REF], it is a nonlinear minimization problem, that may be time consuming to solve. To deal with it, a proposed solution is to replace

u ζ by u * ζ , assuming u ζ u * ζ .
By removing the correction term from the DNN input we can directly deduce the correction term:

∆u ζ = -A -1 ζ δv nd (u ζ ). (22)
Knowing the drone mass, we can get the correction to apply to the thrust reference T ref and to the two angular references φ ref and θ ref to get closer to the linear behavior. The control architecture is summarized in fig. 4. The correction proposed is only done in the position control loop, as attitude control is done considering (3c), which simply corresponds to the rotation of the frame.

C. Stability Analysis

We derive the main lines here to prove that the deep neural network term does not destabilize the quadrotor. Let us consider the following Lyapunov Function:

V (t) = 1 2 η(t) 2 . ( 23
)
Defining the prediction error of the DNN as e η = δ v e -δv nd , and λ = λ min (A -BK η ) be the minimum eigenvalue of the closed loop system matrix.

V (t) = η T (A -BK η )η + η T e η , ≤ λη T η + η • e η , V (t) ≤ 2λV + √ 2V e η . (24)
Assuming e η is bounded by a constant ε η and defining

W = √ V then Ẇ = V 2 √
V it follows that:

Ẇ ≤ λW + ε η √ 2 . ( 25
)
The solution to (25) given W (0) = z0 √ 2 is given by (26):

W (t) ≤ 1 √ 2 ( ε η λ + z 0 )e λ(t-t0) - ε η √ 2λ . ( 26
)
We obtain the following bound:

η(t) ≤ η 0 e λ(t-t0) + ε η λ (e λ(t-t0) -1). ( 27
)
It follows that the solution exponentially converges to a ball, centered at η 0 and with a radius of εη λ . It is directly proportional to the error made by the neural network.

V. SIMULATION AND EXPERIMENTAL RESULTS

In this section, DNN training and validation are explained (Section V-A). Controllers are next tested in ROS/Gazebo simulation environment (Section V-B) and in experimental flights (Section V-D) in our motion capture room.

We work with ROS combined with Gazebo simulation environment to validate our methodology. To perform angular velocity and low level control, PX4 autopilot is used [START_REF]4-autopilot software[END_REF]. Python along with the Keras API [START_REF] Chollet | Keras[END_REF] are used to shape data collected from ROS, named rosbag, to build and to train the DNN. The code used for the simulations and experiments is made open access and available by following the given GitHub link.

A. Training and Validation Database

The proposed Deep Neural Network in Section III is trained with a database created using the linear cascaded control from Section IV-A, for both simulation and experimental tests (fig. 1, step 2 & 3). It is composed of several basic movements as translations in different planes, steps in all directions, rotations and spiral motions (see fig. 5). The idea is to cover all the possible motions of the drone, to better learn the whole dynamic.

The entire state X, the control inputs U and the linear accelerations Ẍ are supposed to be accessible, either by measurements (IMU) or obtained using filters like Extended Kalman Filter (EKF). To validate learned model, one-third of the database is used as validation set. Training and validation sets cover a total of 30 minutes of flight. We propose here to fly with less data as possible to show that we can quickly and correctly learn the error dynamics.

For training process, we used Nesterov Adam optimizer also known as Nadam (Reader is referred to [START_REF] Dozat | Incorporating Nesterov Momentum into Adam[END_REF] for more information about the algorithm). The following hyper-parameters are set for the DNN training: Learning rate l r = 0.05, µ = 0.9, ν = 0.999, batches of size 256 are used to evaluate gradients and a maximal number of 1000 epochs is defined. An epoch refers to the number of times the algorithm has been run through the whole training data-set. In addition, during training process elements are shuffled and the learning rate is decreased by a factor of 2 3 when there is no improvement on the loss function. All those hyper parameters were tuned empirically in order to reach a good compromise between speed and convergence of the algorithm.

We propose to include knowledge of the quadrotor dynamics to the DNN eq. ( 11) by restricting the input variables, using eq. ( 2b): one provides thrust, Euler angles and linear velocities. Sine and cosine of angles are directly provided to ensure that for the DNN 0 and 2π are equivalent, in the same way as in [START_REF] Bansal | Learning quadrotor dynamics using neural network for flight control[END_REF]. Inputs and outputs of the DNN are given by eq. ( 28) and ( 29):

E DN N = z c φ s φ c θ s θ v x v y v z u batt T T , (28)
S DN N = ẍ ÿ z T . ( 29
)
In order to better tackle the ground effect, z component is added to the neural network inputs. The same is applied for battery voltage to handle battery discharge and loss of thrust. Other position components are not included as they do not appear in eq. (2b).

B. Simulations Results

We show results of the cascaded linear controller and it enhanced version with the DNN correction term in ROS/Gazebo simulation environment (fig. 1,step 4). Presented results are given on an unseen trajectory reference, that is to say that they are neither part of the training set nor part of the validation set, see fig. [START_REF] Hunt | Neural networks for control systems -A survey[END_REF]. A third controller, a cascaded linear controller with integral action term is also presented for comparison. In order to be compared correctly with the two other controllers presented, the gains of the controller with integral action are chosen in such a way that the closed loop behavior approaches the ideal behavior without integrator, i.e. the behavior of the LQR controller. Thus the gains proportional, integral and derivative gains are chosen so that the mean squared error between the integral LQR and the non-integral LQR is the lowest. To find these gains a brute force approach was used, requiring several hours of calculations.

When the linear cascaded law is used, reached positions reveal tracking errors and slower time response (blue line). Tracking errors can be explained by the absence of integrator in the initial linear controller. Slower time response is due to the linear control design on a nonlinear model. On tests trajectories, tracking errors with the linear controller are around 7 cm on x and y component and for z component it is up to 25 cm. However, we expected none of those effects if the drone behaves exactly like its linear model (green dashed line).

Input [p q r] [T ] [x y z] Input u [p q r] [T ] Inputs [T ] [u batt ] [vxvy vz ] [z c φ s φ ] [c θ s θ ]
Output ∆u

[δφ δθ δψ] [δT ]
Outputs Those undesired errors, due to nonlinear dynamics of the quadrotor, aerodynamic effects, coupling and unmodelled actuator dynamics (open loop Electronic Speed Controller) are mostly corrected with the proposed corrected DNN Enhanced Linear Cascaded Controller. We can notice that the drone behavior (red line) is closer to the expected linear one: the DNN term has significantly decreased errors. Mean squared errors on this test trajectory are presented in tab. I for each axis. The integral controller tuned to be as close as possible to the linear quadratic controller also corrects the behavior but presents a slower time response.

C. Experimental Setup

In our experiments we use a Holybro Kopis 2 drone. The on-board flight controller is a Kakute F7 where PX4 firmware is installed. Internal PID controller is configured for angular rate loop control. Our drone is connected to a ground station via a wifi bridge that is sending our command signal at 100Hz. Space positions are obtained from the motion capture system, also running at 100Hz. The experimental setup is resumed in fig. 7.

D. Experimental Results

We present results of our experiments using section V-C setup. As for simulation, we first created a database using the linear cascaded controller on the Kopis drone. We collected data with several scenarios using different batteries, including new and old batteries to better cover their discharge profiles. After training our DNN, we tested controllers on an unseen scenario, see fig. 8. Also, to see the impact of battery discharge profile on the controller law, we used a half-full battery for each controller test. The blue line in fig. 8 represents the controller without correction nor integrator: the behavior presents tracking errors. During this flight case, using the linear law, we can notice the huge impact of the battery discharge which is added to the tracking error. At starting, the tracking error is equal to 30% and is increased up to 35% at the end. As expected, the DNN Enhanced Controller completely removed tracking errors and increased flight speed compared to the linear controller. Here the DNN completely compensates the error over the entire test scenario and learned the battery discharge profile as we introduce to the DNN battery discharge tension as an input. Thus, the drone flight behavior is much closer to expected perfect linear dynamic (green dashed line). As for simulation results, the proposed controller is faster and closer to the linear expected behavior than the integral controller. Mean squared errors for experimental test trajectory are presented in tab. II. We invite readers to watch the video of the experiments attached to the paper.

Near Ground Test Comparison

To better see improvements of the DNN Enhanced Linear Cascaded Controller compared to the Cascaded Controller with an integrator, we tested both controllers in a near ground flight case (see fig. 9), here following a lemniscate at 15cm from the ground. In this type of scenario, ground effect strongly affects the drone general behavior. The same deep neural network, as in the previous part, is used as it was trained on many different scenarios including near ground flights. It is able to correct the behavior and get closer to the linear expected behavior in this highly perturbed situation.

Here the integral controller is not rejecting ground effect disturbance, as this effect is not a constant perturbation. The neural term is able to catch those dynamical perturbations and the proposed DNN Enhanced controller, as shown in fig. 10, performs better than the integrator in all three axis.

VI. CONCLUSION AND FUTURE WORK

Rather than tuning a complex nonlinear algorithm, we propose the design of an easy to tune the linear cascaded controller. Then, a learning process is used allowing to correct this basic controller to reach performances initially fixed on a linear quadrotor model. Our experiment shows that with a small amount of flight data collected, the proposed DNN enhancement enables the drone behavior to be closer to a linear one. It allows to reduce undesired overshoots, slow time response and tracking errors. The advantage of this approach is that it is plug-and-play and can be applied to many UAVs present in the industry without the need to find a fine model of the drone. Future work aims at making experiments with online learning and correcting progressively during flight. It can also be extended to learn external perturbations, such as the wind, frame or blades issues, battery unpredicted behavior, etc. Finally, we plan to do a further ablation study to determine the optimal network structure for embedding it in the micro-controller.

Fig. 1 :

 1 Fig. 1: Proposed methodology to fit expected linear behavior on the quadrotor nonlinear dynamics using DNN correction.

Fig. 2 :

 2 Fig. 2: Euler angles references. Inertial frame (left), intermediate frame (middle) and body fixed frame (right).

Fig. 3 :

 3 Fig. 3: Deep Neural Network Architecture: a feed-forward network, with two hidden layers of 64 units, with ReLU function as activation.

  such that we reach the desired closed loop behavior with K ζ . Where gain K ζ is obtained through a linear quadratic regulation synthesis (LQR) and φ θ ψ T ref is the attitude target coming from the position and speed controller. Position Speed controller: Let us define u ζ := [ φ θ T ] T , x 1 := [ x y z ] T and x 2 := [ vx vy vz ] T . Then using (9):

Fig. 4 :

 4 Fig. 4: Proposed Cascaded Controller Architecture for our quadrotor. Attitude and position & speed control loops are running at 100Hz. The DNN correction is learned offline and then is used to correct flight performance with respect to the expected linear dynamic.

Fig. 5 :

 5 Fig. 5: Example of trajectories : step, circles, spirals, etc. included in the database. Distances are in meters.

Fig. 6 :

 6 Fig.6: Evolution of positions for each controller in ROS/Gazebo Simulation Environment test case scenario. The positions of the drone using the cascaded linear controller are represented in blue, with the integrator in yellow, and using the enhanced controller with DNN in red. Green dashed line represents the expected linear behavior on the linear model using cascaded PD linear controller only.

Fig. 7 :

 7 Fig. 7: Experimental setup.

Fig. 8 :

 8 Fig. 8: Evolution of quadrotor positions for each controller, in real flight using our motion capture room (MOCA).

Fig. 9 :Fig. 10 :

 910 Fig. 9: The quadrotor flying close to the ground. Trajectories are displayed in fig. 10

TABLE I :

 I Mean Squared Error (MSE) on each axis for presented scenario (fig.6) for the three controllers in ROS/Gazebo Simulation Environment.

	Computed Command	Motor controller ESC	Secutiry (Failsafe)
	input [ p q r T]			Manual mode commands
			PX4 Flight controller
			Kopis Drone
			Live Odometry (Rostopic)
	Live Position
	Callback
	Vicon		Ground station
	Position captation system	Autonomous Offboard Mode	Compute nominal and corrected command input	Manual Onboard Mode

TABLE II :

 II Mean

Squared Error (MSE) on each axis for presented real flight scenario (fig.

8

) for the three controllers.