Khiem-Van Truong
email: truong@onera.fr

Modeling Aerodynamics Including Dynamic Stall for Comprehensive Analysis of Helicopter Rotors

To fulfill the objective of design tools, comprehensive analysis codes need to be capable of providing both accurate and timeefficient predictions of rotor airloads. The coupling of high-fidelity computational fluid dynamics CFD with comprehensive analysis codes has showed improved predictions of airloads, but such approach is very expensive in terms of CPU time and cannot be used for design. The aerodynamics environment about helicopter rotors is very complex, encompassing subsonic to transonic flow with unsteady, stalled regimes and 3-D effects. Semi-empirical models of dynamic stall were created for modeling unsteady aerodynamics including stalled flow. Most of them, developed from the 1970s through the 1990s, were found unsuccessful to reproduce experimental results on the UH-60A helicopter. Most dynamic stall models also suffered problems of numerical convergence. Thus, there are two levels of difficulties for semi-empirical models: providing good physics description of aerodynamics and ensuring numerical convergence when implemented in comprehensive analysis codes. The present communication is concerned about the revision of the "ONERA Hopf Bifurcation model". The model takes into account various aerodynamic phenomena, unsteady behavior before stall onset, stall delay, vortex-shedding phenomenon and boundary-layer effects that have been overlooked. 3-D effects are not still well investigated and their complexity is accounted for by sweep effects, rotation effects and transonic tip relief effects. For the extreme aerodynamic condition of high-speed flight, a correction outside of the stall model, is suggested for the induced velocity. The implementation of the model into a comprehensive analysis code ensuring the numerical convergence is presented. The experimental results obtained in the wind tunnel S1 of Modane (France) in 1991 on the rotor 7A are considered for illustration of the capabilities of the stall model. There are three test conditions considered: high-speed test point, high-thrust test point with light stall and high-thrust test point with deep stall. Since tests were made in wind tunnel, corrections for test environment involving wind tunnel walls and the test stand, and for Reynolds effects are necessary. The prediction of the airloads and structural loads of all the test points considered is in reasonable agreement with the experimental results and requires very low CPU time demands.

NOTATION

3) in the regime of attached flow C 2,i (i = 1, 2, 3) in the regime of separated flow M Mach number R e Reynolds number α aerodynamic incidence angle α cr critical stall angle r n normalized radial distance R blade radius Λ deg sweep angle, expressed in degrees

INTRODUCTION

To fulfill the objective of design tools, comprehensive analysis codes need to be capable of providing both accurate and time-efficient predictions of rotor airloads. The coupling of high-fidelity computational fluid dynamics CFD with comprehensive analysis codes has showed improved predictions of airloads on the UH-60A [START_REF] Potsdam | Rotor Airloads Prediction Using Loose Aerodynamic/Structural Coupling[END_REF], but such approach is very expensive in terms of CPU time and cannot be used for design. The aerodynamics environment about helicopter rotors is very complex, encompassing subsonic to transonic flow with unsteady, stalled regimes and 3-D effects. Semi-empirical models of dynamic stall were created for modeling unsteady aerodynamics including stalled flow. Most of them were developed from the 1970s through the 1990s ( [START_REF] Tran | Semi-empirical model for the dynamic stall of airfoils in view of the application to the calculation of responses of a helicopter blade in forward flight[END_REF], [START_REF] Petot | Differential Equation Modeling of Dynamic Stall[END_REF], [START_REF] Leishman | A semi-empirical model for dynamic stall[END_REF], [START_REF] Truong | A 2-D dynamic stall model based on Hopf bifurcation[END_REF]) but failed to provide good predictions of experimental results on the case of the UH-60A in stalled regime [START_REF] Nguyen | Evaluation of Dynamic Stall Models with UH60A Airloads Flight Test Data[END_REF]. Applied to other stalled test [START_REF] Yeo | Calculation of Rotor Performance and Loads Under Stalled Conditions[END_REF], most dynamic stall models also suffered problems of numerical convergence. Thus, there are two levels of difficulties for semi-empirical models: providing good physics description of aerodynamics and ensuring numerical convergence when implemented in comprehensive analysis codes. The present communication is concerned about the revision of the "ONERA BH model" [START_REF] Truong | A 2-D dynamic stall model based on Hopf bifurcation[END_REF], to be named from now on "ONERA Hopf bifurcation model". This model received initially a positive review [START_REF] Snel | Review of the Present Status of Rotor Aerodynamics[END_REF], but was later severely criticized for its account of shedding-vortex phenomenon in the airfoil wake by a Van-der-Pol type equation. The experimental results obtained in the wind tunnel S1 of Modane (France) in 1991 on the rotor 7A are considered for the illustration of the capabilities of the revised stall model. There are three flight conditions considered: the high-speed test point 312, the high-thrust test point 293 with light stall and the high-thrust test point 596 with deep stall. The computations are done with the comprehensive analysis code ROTOR of ONERA equipped with the revised dynamic stall model. Many studies have been made for the highspeed test point ( [START_REF] Pahlke | Chimera Simulations of Multibladed Rotors in High-Speed Forward Flight with Weak Fluid-Structure-Coupling[END_REF], [START_REF] Beaumier | Weak and strong coupling between the elsA CFD solver and the HOST helicopter[END_REF]), but no correlation studies have been published for the second experimental test point, until recently ( [START_REF] Ortun | Rotor Loads Prediction on the ONERA 7A Rotor using Loose Fluid/Structure Coupling[END_REF], [START_REF] Surrey | Advanced modeling of the 7A structural rotor blade dynamics[END_REF]). No study has been published for the third data test point. In the next chapter, the modeling of aerodynamics is detailed, first for dynamic stall phenomenon. Aerodynamic environment for the rotor blade is inherently 3-D; so, it is necessary to introduce various effects, sweep effects, rotation effects and transonic tip relief effects. For the extreme aerodynamic condition of high-speed test point, a correction of the induced velocity, outside of the stall model, is proposed. Corrections for the test environment involving the wind tunnel walls and the test stand and for Reynolds effects are also necessary. In the third chapter, the implementation of the dynamic stall model in the comprehensive analysis code of ONERA is described. In the fourth chapter, the application of the comprehensive analysis code equipped with the dynamic stall model is made to the three different test points of the rotor 7A. Concluding remarks will be drawn in the last chapter.

MODELING AERODYNAMICS

Aerodynamics about the rotor results from the interaction of the incoming flow with the rotor and can be viewed as composed of a global flow that extends from far upstream of the rotor disc to far downstream, and a local flow about the rotor blades. This separation is used in comprehensive analysis to solve the aerodynamics using different methods for the global flow field and for the local flow field. Nevertheless, the two flow regions are interconnected and this feature is taken into account, the inflow conditions on the rotor blades are determined by the global flow and the aerodynamic forces acting on the blades influence the global flow. The global flow is evaluated by momentum or wake methods. The aerodynamic forces on the blades can be taken from wind tunnel experiments (generally completed by computation) and presented as tables of non-dimensional lift, drag and pitching moment coefficients with respect to aerodynamic incidence angles and Mach numbers. A number of corrections are applied to the 2-D values through models of dynamic stall. The publicized objective of dynamic stall models has been the modeling of dynamic stall aerodynamics. Actually, it involves also modeling of unsteady aerodynamics prior to the occurrence of dynamic stall. This study will be concerned with various improvements of the "ONERA Hopf Bifurcation" stall model.

Modeling Dynamic stall phenomenon

In the limited extent of this communication, it will not be possible to detail all the aspects of the dynamic stall model.

Let us briefly summarize that the flow is considered as attached or still attached (due to stall delay) and separated flow. The aerodynamic coefficients C 1,i (i = 1, 2, 3 for lift, pitching moment and drag respectively) can be separated into attached and separated components:

C i = C 1,i +C 2,i (1)
In the first regime, the aerodynamic coefficients can be separated into impulsive and circulatory components:

C 1,i = C impulsive 1,i +C circulatory 1,i (2)
The values of C 1,i vary around the static values. The latter values are obtained from the fitting of experimental polars according to analytical relations similar to those of the Leishman-Beddoes model ( [START_REF] Beddoes | Representation of Airfoil Behavior[END_REF], [START_REF] Leishman | A semi-empirical model for dynamic stall[END_REF], [START_REF] Truong | A 2-D dynamic stall model based on Hopf bifurcation[END_REF]).

The present dynamic stall model distinguishes from all the other models by its analysis of the separated flow. To illustrate the reasoning, consider Fig. 1 that shows the behavior of the lift coefficient during a periodic cycle of the pitching angle of an airfoil. For the state between the points 1 and 2, the flow remains attached. It is still attached up to the point 3 due to flow delay [START_REF] Ericsson | Fluid mechanics of dynamic stall Part I. Unsteady flow concepts[END_REF], although the pitching angle has increased beyond the value of the static critical angle of stall. From the point 3 to point 4, there is the stall development stage that is usually defined as uniquely the dynamic stall vortex formation. Actually, there are two flow phases ( [START_REF] Jumper | Toward Unsteady Lift Augmentation: An Assessment of the Role of Competing Phenomena in Dynamic Stall[END_REF], [START_REF] Mulleners | Dynamic stall development[END_REF]). In the first phase, there is formation of small scale co-rotating vortices that emerge in the shear layer, according to Mulleners et al. [START_REF] Mulleners | Dynamic stall development[END_REF]. Then the shear layer goes through a primary instability stage, this instability is followed by a vortex formation stage, characterized by the roll-up of the shear layer into the large-scale dynamic stall DS vortex. Following its formation, the DS vortex convects towards the trailing edge of the airfoil. Between point 4 and point 7, the flow remains separated and there is formation of multiple vortices, leading to various maximum of the lift coefficient. In Fig. 1, there are only two secondary vortex, leading to two peak maximums at point 5 and 6 respectively. From point 7, the flow begins its reattachment regime.

According to the present stall model, the vortex-shedding phenomenon is simulated as a self-excited oscillator of Vander-Pol type equation when the aerodynamic incidence angle (or angle of attack) keeps growing above the dynamic critical stall angle:

d 2 C 2,i dt 2 -ω S (β + i -γ + i C 2 2,i ) dC 2,i dt + ω 2 S C 2,i = -E + i ω S d|α| dt (3
) and by a damped oscillator when the flow is reattached:

d 2 C 2,i dt 2 -ω S β - i dC 2,i dt + ω 2 S C 2,i = 0 ( 4
)
where ω S , β + i , γ + i , E + i and β - i are constants. The above equations have been simplified from the original modeling [START_REF] Truong | A 2-D dynamic stall model based on Hopf bifurcation[END_REF], there is only one term proportional to d|α| dt in the second member of Eq. 3. The Van-der-Pol equation has been used for vortex-shedding phenomenon past bluff bodies [START_REF] Blevins | Flow-Induced Vibration[END_REF]. The modeling of Eq. [START_REF] Blevins | Flow-Induced Vibration[END_REF][START_REF] Bousman | Putting the Aero Back Into Aeroelasticity[END_REF] does not raise interest to other researchers, except a few ones [START_REF] Snel | Review of the Present Status of Rotor Aerodynamics[END_REF], due to its non-linear behavior that leads to numerical difficulties for comprehensive analysis codes. However, according to the physical scenario described above, the excess of unsteady aerodynamic coefficient associated to the primary instability of the shear layer needs to be accounted for. Its contribution is in the same order as from the vortex formation. In the model published previously [START_REF] Truong | A 2-D dynamic stall model based on Hopf bifurcation[END_REF], this effect is not clearly stated but it is taken into account by the term ∂C 1,i /∂α (α: aerodynamic incidence angle), obtained from the response indicial approach. Such term appears also in the ONERA-EDLIN model (first term in the right hand side of Eq. 7 of [START_REF] Petot | Differential Equation Modeling of Dynamic Stall[END_REF]). The formulation in terms of the derivatives ∂C 1,i /∂α requires precise knowledge of the aerodynamic coefficients in the angle region around the critical stall angle that is not usually available. So, it is proposed to delete the derivative component in the first equation governing the aerodynamic coefficients and instead add the following equation:

d CBL dt = λ BL (C equil BL -CBL ) (5)
where:

C equil BL =    0.0 : i f non -stall -c 0 BL • d f Bed dt • c/2 v n : otherwise (6)
where c 0 BL is a constant, f Bed denotes the separation function introduced by Beddoes [START_REF] Leishman | A semi-empirical model for dynamic stall[END_REF], v n is the fluid velocity normal to the blade section and c is the chord.

3-D effects

The aerodynamic environment about the rotor blades is inherently 3-D. The knowledge of 3-D effects is still partial, there are ongoing studies based on experiments and CFD techniques. In comprehensive analysis, there are already semi-empirical corrections of 3-D effects. Harris [START_REF] Harris | Rotor High Speed Performance, Theory vs Test[END_REF] is the first to propose corrections for radial flow based on sweep angles. This correction leads basically to the delay of the critical stall angle. This type of correction improves the predictions, but it is not known that the variable in play is the instantaneous sweep angle or some delayed value. In wind turbines, it is observed that the normal force is increased in inboard blade sections r n < 0.5 (r n : normalized radial distance). Such effects are called "rotation effects" and appear to be due to centrifugal and Coriolis forces [START_REF] Breton | A Study on Rotational Effects and Different Stall Delay Models Using a Prescribed Wake Vortex Scheme and NREL Phase VI Experiment Data[END_REF]. Various corrections have been proposed for the aerodynamic coefficients [START_REF] Breton | A Study on Rotational Effects and Different Stall Delay Models Using a Prescribed Wake Vortex Scheme and NREL Phase VI Experiment Data[END_REF], but none provides good correlation for the predictions with the experiments. In the present model, the following correction is proposed and it is limited to the normal force:

C rot n = 0.0 : i f α ≤ α 3-D
c 0 rot • sin(α -α 3-D cr ) • f 0 (r n ) : otherwise (7
) where f 0 (r n ) denotes a function decaying rapidly to zero as r n increases from 0.5:

f 0 (r n ) = 1.0 + tanh((0.5 -r n )/0.10)) (8)
There is also a new correction introduced in the stall model, associated with finite span blade effects. The phenomenon is described in Sitaraman thesis [START_REF] Sitaraman | CFD based unsteady aerodynamic modeling for rotor aeroelastic analysis[END_REF] as transonic tip relief effects appearing at the blade tip, leading to an important change of the pitching moment coefficient over a blade region of about 30% radius span. It is proposed to simulate the behavior of the pitching moment by the following equations:

C 3-D m = C 2-D m • f 1 (r n ) +C Tip m • f 2 (r n ) (9)
where f 1 (r n ) denotes a function decaying rapidly at the blade tip and f 2 (r n ) decays towards inboard blade sections:

f 1 (r n ) = tanh((1.0 -r n )/0.10)) f 2 (r n ) = 0.50 • (1.0 + tanh((r n -0.82)/0.05)) (10)
C Tip m = C T 1 m +C T 2 m C T 1 m =(c 0 m +c 1 m •α+c 2 m •α• α +c 3 m •α 3 ) •g(M)•c 4 m C T 2 m = g(M) • h(Λ deg ) • c 5 m      (11)
where α denotes the aerodynamic incidence angle, M the Mach number, Λ deg the sweep angle in degrees,

c 0 m , c 1 m , c 2 m , c 3 m , c 4
m and c 5 m are constants. The functions g and h are defined by the following relations:

g(M) = 0.5 + 0.25 • (1.0 + tanh((M -0.7)/0.1)) h(Λ deg ) = 0.5 • (1.0 -tanh(Λ deg /20.0)) (12)
According to the analytical expression for g(M), the effects increase with Mach numbers but remain limited when M exceeds 0.7. The transonic tip relief effects depend also on the sweep angle, the dependence follows the law in "tanh" that limits the effects when the sweep angle exceeds 20 degrees. The variation C T 1 m with the aerodynamic incidence angle is suggested by the values computed by Sitaraman. The additional component C T 2 m that depends on sweep angles is needed to create the asymmetry, with respect to azimuth, of the values of the normal force about its minimum value at around the azimuth of 100 degrees. The behavior of the pitching moment, described by Eq. (9-12), has been proposed, guided by physical considerations. It is possible to verify the suggested behavior by studying swept blades with high-fidelity numerical tools.

Aerodynamic correction for the high-speed flight

It is well known that comprehensive analysis fails to predict the airloads of rotors in high-speed flights. Bousman [START_REF] Bousman | Putting the Aero Back Into Aeroelasticity[END_REF] raised the two main deficiencies of comprehensive analysis consisting in the predictions of the azimuthal position of the negative lift force in the advancing side of the UH-60A and the under-prediction of the pitching moment, leading to the under-prediction of the pitch link loads. These two issues have been resolved with the coupling of comprehensive analysis with CFD techniques [START_REF] Potsdam | Rotor Airloads Prediction Using Loose Aerodynamic/Structural Coupling[END_REF]. The same issues occur for the rotor 7A and they have been recently resolved ( [START_REF] Ortun | Rotor Loads Prediction on the ONERA 7A Rotor using Loose Fluid/Structure Coupling[END_REF], [START_REF] Surrey | Advanced modeling of the 7A structural rotor blade dynamics[END_REF]). The principal difference between computation made with comprehensive analysis standalone and comprehensive analysis coupled with CFD techniques resides in the consideration of viscous effects by advanced aerodynamics methods and ignored in low-order methods.

The principal feature resulting from viscosity consideration is the phase shift in azimuthal position of the normal force for the blade sections in the tip region. To improve the capabilities of global inflow models, the following correction to the induced velocity is proposed:

v corr ind = c 0 ind • Ω • R • (-r n ) • cos(2.0 • ψ + φ 0 ) (13)
where c 0 ind and φ 0 are constants, Ω denotes the rotation speed, ψ the azimuth and R the blade radius. The effects of such correction on typical values of the induced velocity at the blade tip can be visualized in Fig. 2. The maximum of the induced velocity is shifted from the azimuth at 90 degrees to higher values. Since the angle of aerodynamic incidence depends on the induced velocity, its minimum is also shifted and therefore the position of the minimum of the normal force to higher value of azimuth.

Wind tunnel effects

Since the 7A rotor is placed in the wind tunnel S1 of ONERA, the aerodynamic environment about the rotor blades is perturbed by the presence of the hub and by the wind tunnel walls. There is a correction accounting for the rotor hub and the wind tunnel wall effects made by Aerodynamics Department of ONERA in 1993 (Internal report of A. Masson and P. Beaumier), based on a surface panel method. The computation resulted into a correction of the induced velocity on the nodes positioned on the blade sections. Fig. 3 shows the values of the induced velocity normalized to the incoming flow velocity (according to the direction of axis x) for the shaft angle α q = -8.0 degrees. Only a half azimuthal region is drawn, the other half can be deduced by symmetry. It can be seen that the incoming flow is deviated by the rotor hub. Another important effect introduced by wind tunnels concerns Reynolds numbers. Indeed, the polars (look-up tables of aerodynamic coefficients) are provided for high Reynolds numbers, in the case of the rotor 7A for the ratio R e /M around 8 millions. The value of Reynolds numbers affects not only the critical stall angle but also the dependence of the values of the aerodynamic coefficients versus the aerodynamic incidence angles. Typically, the function f α of the lift force versus the angle of attack is concave at high Reynolds numbers and it becomes more rounded as the Reynolds number decreases [START_REF] Jacobs | Airfoil section characteristics as affected by variations of the Reynolds number[END_REF]. As there is no available data for the airfoils OA209 and OA213 that composed the 7A blade, the data on the family of NACA airfoils is used [START_REF] Jacobs | Airfoil section characteristics as affected by variations of the Reynolds number[END_REF]. The correction is limited to only the values of critical stall angle and not on the shape dependence f α . Neglecting the variation of the normal force versus the angle of attack with Reynolds number will imply the assumption of keeping the same shape dependence, this will lead to lower values of lift force at high angles of attack (greater than the critical stall angle). The experimental values of the critical stall angles for the airfoils NACA0009, NACA0012, NACA0015 and NACA2412 are reported in Fig. 4. These experimental values are fitted with the analytical expression of the following form:

α cr = c 0 R e + c 1 R e • tanh((R e • 1.0E -6 -c 2 R e )/c 3 R e ) (14)
where c 0 R e , c 1 R e , c 2 R e and c 3 R e are constants.

IMPLEMENTATION OF THE STALL MODEL

The Hopf bifurcation model is implemented as a module coupled to the comprehensive analysis code ROTOR of ONERA. The version of ROTOR used is a recently developed code, with predicting capabilities for the whole helicopter, instead of an isolated rotor, as in the old version developed in the 1990s. The coupling is similar to the loose coupling between comprehensive analysis code and CFD code, but with different variables for the communication between the two codes. The comprehensive analysis code sends to the stall module the values of the aerodynamic conditions on the blade nodes during a rotation period, i.e. the aerodynamic incidence angle, the pitch angle and the fluid velocity adjacent to the blade. Based on these values, the stall module will calculate the values of the airloads (normal force, pitching moment, drag and radial forces) on the blade nodes for the whole rotation period. It iterates as many times as needed to attain a periodic solution. Since the aerodynamic forces are damped rapidly, it is sufficient to run the computation for three periods. However, for the particular study of flow delay with slow time delay, it was necessary to develop a method based on Poincaré maps for calculating the periodic solution [START_REF] Parker | Practical numerical algorithms for chaotic systems[END_REF]. The method relies on the fact that the periodic solution of the mechanical system, if it exists, will be a fixed point for the intersection of the trajectory of the mechanical system with the phase plane at times equal to multiple integers of the period. Accordingly, the problem is reduced to the search of a fixed point by Newton's method. The module is built for allowing the study of separate behaviors of the aerodynamic flow: quasi-steady "s", unsteady without stall "u", delayed flow "d", flow with boundary-layer effects "b" and flow with vortex-shedding phenomenon "v". It permits any combination of the various behaviors of the flow cited above. Such capabilities of the code permits separately validation of the parameters associated with each particular flow behavior. The general case with stall is determined buy the option "udbv". As stated in Eq. 2, the aerodynamic coefficients are decomposed into impulsive and circulatory components. The equations governing the unsteady behavior of the impulsive component are taken from [START_REF] Johnson | Rotorcraft Dynamics Models for a Comprehensive Analysis[END_REF]. based on thin airfoil theory. According to Johnson [START_REF] Johnson | Rotorcraft Dynamics Models for a Comprehensive Analysis[END_REF], the circulatory component is taken into account by the vortex wake model. The parameters for the delayed flow are provided by Beddoes [START_REF] Beddoes | Representation of Airfoil Behavior[END_REF]. The remaining unknown parameters are those characterizing the boundary-layer effects and the vortex-shedding phenomenon.

To further simplify the model, only the equation governing the pitching moment for the behaviors of boundary layer and for vortex shedding are solved, the components of the normal and drag forces are obtained by a multiplicative constant (4× for the normal force, and 1.6× for the drag force). This choice is made from the analysis of the experiments on the NACA0012 airfoil. It could be changed accordingly to the airfoil studied.

For any choice of the behavior of the aerodynamic flow, the aerodynamic forces on the blade sections of the rotor are governed by a system of ordinary differential equations ODE.

With the system of ODE involving the stall phenomena (Eq. 3-4), the flow equations will contain two different regimes, the first involving the unsteady flow without separation that varies slowly upon time while the regime involving the vortex-shedding phenomena has a rapid time-varying behavior. The system of ODE with two different behaviors with respect to time are qualified as "stiff". There are two types of ODE solvers,"non-stiff" and "stiff". The procedure for solving the system of ODE of the Hopf bifurcation model would be the use of "non-stiff" solver before stall onset and "stiff" solver at stall onset. Fortunately, there exists a solver specially developed for such circumstances, "LSODA" by Hindmarsh and Petzold in the 1980s ( [START_REF] Hindmarsh | A Systematized Collection of ODE Solvers[END_REF], [START_REF] Petzold | Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations[END_REF]), switching automatically from non-stiff to stiff solver whenever it is needed. The use of a unique solver "stiff" or "non-stiff" is not time-efficient and may lead to numerical divergence.

The system of equations governing the mechanical behavior of the helicopter is non-linear with the addition of the equations governing the aerodynamic forces on the rotors (principal and tail). The search for an equilibrium solution is not an easy numerical task. The comprehensive analysis code ROTOR uses the solver hybrd of the package MINPACK [START_REF] Moré | User Guide for MINPACK-1[END_REF]. It is based on a modification of Powell's hybrid algorithm. Although developed in the 1980s, it is quite robust and has solving capabilities comparable with other solvers. Its principal weakness consists in the request of the initial solution not too-far from the final solution; otherwise, it may diverge or take many iterations to converge.

APPLICATION TO THE ROTOR 7A

In 1991, ONERA made a series of experiments in the wind tunnel S1 at Modane on the rotor 7A to provide test cases for analysis. This rotor is articulated and composed of four blades with rectangular planform and a radius of 2.1 meters. The experiments were made in various configurations, from attached flow regime to stalled flow regime.This study will consider three test conditions: high-speed test point 312, high-thrust test point 293 with light stall and high-thrust test point 596 with deep stall. The last configuration was obtained by slowing the rotor speed, from its nominal speed of about 1020 rpm to about 700 rpm. Many attempts to improve analysis were made for getting a good correlation between predictions and experiments ( [START_REF] Pahlke | Chimera Simulations of Multibladed Rotors in High-Speed Forward Flight with Weak Fluid-Structure-Coupling[END_REF], [START_REF] Beaumier | Weak and strong coupling between the elsA CFD solver and the HOST helicopter[END_REF]) and only recently [START_REF] Ortun | Rotor Loads Prediction on the ONERA 7A Rotor using Loose Fluid/Structure Coupling[END_REF], the coupling of comprehensive analysis with CFD techniques meet success for the test points 312 and 293. The deep stall test point 596 is not yet considered for coupled analysis.

The computation is made for all the test points with the global flow simulated by the prescribed wake model, uniform inflow or dynamic inflow models are not used since they do not provide good values of angle of attack in the tip region of the blade.

High-speed test point 312

The high-speed test point 312, characterized by the advance ratio µ = 0.4 and the lift force Z = 12.0 was the case of choice for testing the capabilities of analysis in France and Germany. The comprehensive analysis code showed the same deficiencies in predictions, as for the UH-60A rotor [START_REF] Bousman | Putting the Aero Back Into Aeroelasticity[END_REF], i.e. the shift in azimuthal position of the minimum of the normal force in the advancing side, and the absence of positive values of the pitching moment in the first azimuthal quadrant for the blade sections in the tip region. These two deficiencies are shown in Fig. 5 depicting the normal force and Fig. 6 depicting the pitching moment at various blade sections, for the analysis made under the assumption of quasi-steady aerodynamics. The correction due to the wind tunnel effects is not sufficient to improve the predictions: there is mainly a bump created for the normal force at the azimuth around 80 o for the blade sections of 97%R and 91%R. The remedy to the second deficiency of the prediction is to take account of the transonic tip relief effects on the pitching moment. As explained above, the solution has been to describe the pitching moment according to Eq. [START_REF] Hindmarsh | A Systematized Collection of ODE Solvers[END_REF][START_REF] Jacobs | Airfoil section characteristics as affected by variations of the Reynolds number[END_REF][START_REF] Johnson | Rotorcraft Dynamics Models for a Comprehensive Analysis[END_REF][START_REF] Jumper | Toward Unsteady Lift Augmentation: An Assessment of the Role of Competing Phenomena in Dynamic Stall[END_REF]. The values of the parameters chosen for the point 312 are the following:

c 0 m = -0 • 08, c 1 m = 1 • 081E -2, c 2 m = -1 • 213E -4 c 3 m = 2 • 855E -5, c 4 m = 0 • 12, c 5 m = 0 • 09 (15)
To create the shift of the azimuthal position of the minimum of the normal force, the correction of the induced velocity according to Eq.( 13) is used with the following values of the parameters:

c 0 ind = 0 • 02, φ 0 = 38.0 × (π/180.) (16)
The correction of the induced velocity assumed according to Eq.13 is not large, as seen in Fig. 7, its main effect is to shift the position of the maximum (according to the axis convention of ROTOR) of the induced velocity at around the azimuth of 90 o to a higher value. With the conjunction of wind tunnel corrections, tip effects, rotation effects, unsteady behavior and delayed flow, the predictions of the normal force and the pitching moment are improved. They are not as good as from coupled analysis of comprehensive analysis and CFD techniques but still present much improvement compared to the initial predictions under quasi-steady aerodynamics and at the small expense of about 15 minutes of CPU time on a PC, that expense is very small compared to coupled analysis.

A remark to make is about the very negative values of the pitching moment at inboard blade sections, they are due to the small number of pressure sensors on the blade (as shown in [START_REF] Ortun | Rotor Loads Prediction on the ONERA 7A Rotor using Loose Fluid/Structure Coupling[END_REF]). The improvement brought by to the stall model can be also visualized on the vibratory loads (harmonics 3-8 of the normal force), the stall model provides a better phase to the predictions for all the outboard sections up to the section 50%R where the predictions are degraded (see Fig. 8).

The improvements of airloads lead to improved predictions of the structural loads of flap bending moment and torsion moment. For the latter, the stall model provides the phase shift, although not perfect. The amplitude predicted of the torsion moment is about 20% higher that experimental results. For the chord bending moment, the predictions are not only higher than experimental results but also contain four peaks instead of three as in experiments. It appears that the code ROTOR needs revision for the evaluation of the chord bending moment.

The rotation effects are assumed to provide improved predictions of the normal force at inboard blade section, near the rotor hub. Unfortunately, measurements on the rotor 7A were limited to inboard section at the radius of 0.50%R. Nonetheless, it is possible to estimate indirectly the corrections due to rotation effects. Indeed, to get the trim on the global lift of the rotor, the normal forces have to be equilibrated on all the blade sections. Since rotation effects provide more amplitude for the normal forces at inboard sections, the normal forces at outboard sections are submitted to a decrease in order to maintain the equilibrium; particularly at the section 97%R, the normal force at the azimuth of around 100 o decreases to attain more negative values.

The rotor control angles are the three components of the pitch angle -collective θ, longitudinal cyclic θ 1s and lateral cyclic θ 1c angles -and the shaft angle α q . Calculations have been made under trim objectives of the experimental values of the rotor lift force Z, the rotor propulsive force and Modane flapping laws (zero lateral flapping and longitudinal cyclic equal to minus longitudinal flapping angle). Table 1 lists the values of the rotor control angles from experiments, and from calculations made with quasi-steady aerodynamics, or with wind tunnel corrections or with the stall model. The values of the shaft angle calculated with all the options are reasonably correct. The values of the longitudinal cyclic calculated are different from the experiments, the difference is between 1.0 o and 2.0 o , according to the option of aerodynamics. The worst results are obtained for the lateral cyclic, with about 2.2 o of difference for the option of stall model. The predictions of power P are overestimated of about 8% from the test values.

The old version of ROTOR predicts control angles that are not far from experiments, it is necessary to improve the new version for control angles calculation.

High-thrust test point 293

This test point is characterized by the advance ratio µ = 0.3 and the lift force Z = 19.0 The predictions of this test flight have been recently improved by coupled analysis of comprehensive analysis with CFD techniques, by teams at US. Army and ONERA [START_REF] Ortun | Rotor Loads Prediction on the ONERA 7A Rotor using Loose Fluid/Structure Coupling[END_REF] . However, no computation based on comprehensive analysis standalone with reasonable correlation between predictions and experimental results has been published. Computation made so far fails to detect stall. The cause is due to Reynolds effects, responsible for decreasing the critical stall angle, as shown previously in Fig. 4. With the account of critical stall angles based on the family of NACA airfoils, it is possible to correct the critical stall angle for each blade section on the blade 7A. Through this correction, the comprehensive analysis predicts, with the stall model implemented, the appearance of stall phenomenon for this test point. The parameters of the stall model for this test point are the same as for the previous test point, except:

c 4 m = -0.12 • (0.3/0.4), c 5 m = 0.09 • (0.3/0.4) c 0 ind = 0.02 • (0.3/0.4), φ 0 = 38.0 • (0.3/0.4) (17)
The above constants are taken equal to those of the previous test point with a multiplicative factor proportional to the advance ratio.

The values of parameters related to vortex-shedding phenomenon and boundary-layer effects are chosen as:

ω S = 0.075 • 2.0 • π; β + i = 0.008; γ + i = 1.700 E + i = 0.15; β - i = -3.0; λ BL = 0.2 (18)
The predicted values of the pitching moment in Fig. 11 show negative peak values for all the blade sections, except at the blade section at 97%R. Computation has been made also with another comprehensive analysis code, HOST, developed by Airbus Helicopters, but with no stall model. Its results do not show any stall, due to the fact explained previously. Note that the predicted stall azimuthal region is larger than in experiments. For the normal force shown in Fig. 10, the minimum at around the azimuth of 270 o is captured except for the blade section at 97%R. There is also a singular behavior of the normal force at the azimuth of around 70 o , it has a peak at this value and the peak is very pronounced at this blade section. The amplitude is much greater if the slope of the lift force is chosen to be given by the relationship of Prandlt-Glauert (equal to 1.0/ (1.0 -M 2 ). In the result shown in Fig. 10, its value is chosen not to follow the above relationship but to decrease as M attains 0. Usually, the comprehensive analysis code does not succeed to obtain numerical convergence at such high values of the lift force and calculation is made under fixed control settings, given by experimental values. For both calculation options, the comprehensive analysis code ROTOR equipped with the Hopf bifurcation model succeeds to provide converged numerical solution. For airloads, the calculations provide a good correlation with experiments (Fig. 14). The number of secondary vortex shed over the blade is not correctly predicted for the blade section 82.0%R, but the amplitude of the pitching moment at stall onset is reasonably predicted. The experimental results of the normal force in the fourth azimuthal quadrant are not well reproduced by the stall model, leading to the same deficiencies for the vibratory loads (Fig. 14). The calculations under trimmed conditions provide better predictions for the flap bending moment (Fig. 15) but deviate much from experiments for the outboard sections 55%R and 75%R. They have similar waveform as in experiments, but have less oscillations. For the inboard sections at 30%R and 40%R, the correlation is good except in the range of azimuth [0. o ,60. o ]. For this test with great influence of dynamic stall, the amplitude of the chord bending moment predicted (Fig. 16) does not exceed too much the experimental values. The number of oscillations is visibly the same as in experiments. The predicted values of the torsion moment follow the waveform of the experimental results (Fig. 17). However, the number of oscillations is greater than observed experimentally. This feature will need more analysis. The power prediction is very large compared to experiments, of about 40% of overestimation. One probable cause is the dependence of the lift coefficient versus the aerodynamic incidence angle, it has a shape more rounded when the Reynolds number decreases. As consequences, the lift force should get higher values for the same high value of incidence angle at low Reynolds number than at high Reynolds number.

The trim objective of lift force will then be obtained at lower pitch angle, this equilibrium will lead to lower drag values and therefore lower power.

The importance of stall phenomenon can be seized through the examination of airloads and structural loads of the three test points. As the 7A rotor enters deeper in the stall regime, the pitching moment is observed to rise rapidly with an increasing amplitude at the beginning of the fourth azimuthal quadrant. This variation is the key component influencing the structural loads, particularly the torsion moment behavior that follows more closely the pitching moment (see Fig. 9, 13 and 17). The flap and chordwise bending moments acquire a richer harmonic component (see Fig. 9, 13, 15 and 16).

SUMMARY AND CONCLUSIONS

To attain reasonable agreement of the predictions of the comprehensive analysis with experiments on the rotor 7A, various issues have been solved at two different levels, physics modeling of aerodynamics and code implementation for ensuring numerical convergence: 1. Action of modeling aerodynamics is made at the level of global flow, outside of the stall model, for the account of the extreme aerodynamic condition of the high-speed flight. The proposed correction could be checked independently by CFD methods.

2. The Hopf bifurcation stall model has been revised to emphasize the effects of boundary-layer accompanying the formation of the dynamic stall vortex. This phenomenon is responsible for about half of the variation of the pitching moment at stall onset. This stall model is the only one to treat vortex-shedding phenomenon as a non-linear effect governed by ordinary differential equations.

3. Comprehensive analysis use look-up tables for values of aerodynamic coefficients and thus are time-efficient means for predictions of the aeromechanical behavior of the rotor. These tables provide 2-D values of aerodynamic coefficients. Various corrections are necessary for taking into account the 3-D effects, they are based on sweep effects, rotation effects and transonic tip relief effects.

4. The procedure for the implementation of the Hopf bifurcation model to the comprehensive analysis code ROTOR of ONERA has been presented. The numerical code established is robust and allows convergence in the deep-stall test point of the rotor 7A.

5. The application of the comprehensive analysis code ROTOR equipped with the stall model has been illustrated on three test points of the rotor 7A carried in the wind tunnel S1 of ONERA, involving high-speed and high-thrust test points. For all test points, the predictions of airloads are in reasonable agreement with experiments and request about 35 minutes on a PC for the stalled test point and less for the high-speed test point. Improvement of airloads predictions leads to improved predictions of structural loads.

ACKNOWLEDGMENTS

The author would like to express his thanks to Mr. D.Petot for his explanations about the use and the architecture of the code ROTOR.

  cr and α ≤ 60.0 * (π/180.0)

The 2 -

 2 D values of the pitching moment coefficient C 2-D m decreases rapidly to zero near the blade tip, while the values C Tip m extends to the inboard region of the blade over a length around 0.18 r n . The following equations govern the behavior of C Tip m :

Fig. 1 :

 1 Fig. 1: Lift force of an airfoil pitching in deep stall regime

Fig. 2 :Fig. 3 :Fig. 4 :Fig. 5 :Fig. 6 :Fig. 7 :Fig. 8 :Fig. 9 :Fig. 10 :Fig. 11 :Fig. 12 :Fig. 13 :Fig. 14 :Fig. 15 :Fig. 16 :Fig. 17 :

 234567891011121314151617 Fig. 2: Correction of Induced velocity for high-speed flights

Table 1 :

 1 High-speed flight 312: control angles and power

		θ	θ 1c	θ 1s	α q	P (kW)
	Experiments 10.41 3.43 -3.70 -13.75	88.0
	QS aero.	11.28 0.99 -1.33 -13.48	98.17
	W. Tun. Cor. 12.11 2.16 -1.65 -13.59	92.32
	Stall model	12.11 1.20 -1.82 -13.92	95.10

Table 2 :

 2 High-thrust flight 293: control angles and power

		θ	θ 1c	θ 1s	α q	P (kW)
	Experiments	8.40 3.16 -3.51 -6.70	74.80
	w/o St. Mod. 10.16 1.90 -1.98 -7.90	78.13
	w Stall model 7.74 1.09 -2.20 -7.37	99.92

Table 3 :

 3 High-thrust flight 596: control angles and power under fixed and trimmed control settings (C.T.) As for the high-speed test point, the control angles given by calculations, shown in Table2, do not have a good correlation with experiments. The value of the collective pitch angle calculated by HOST is higher than experiments of about 1.8 o , it is understandable because the aerodynamic model in this code does not provide enough lift force. The ROTOR code provides reasonable values for the shaft angle and the collective pitch angle, but as for the previous test, it gives low value of longitudinal cyclic compared with experiments, and especially very low value of lateral cyclic. For the power prediction, again this code predicts a too-high value of about 32%, compared to experiments.

	8. Coupled analysis


Rafael Granero-Belinch Ón

GLOBAL SOLUTIONS FOR A HYPERBOLIC-PARABOLIC SYSTEM OF CHEMOTAXIS

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

Introduction

In this note we study the following system of partial differential equations

∂ t u = -Λ α u + ∇ • (uq), for x ∈ T d , t ≥ 0, (1)
∂ t q = ∇f (u), for x ∈ T d , t ≥ 0, ( 2
)
where u is a non-negative scalar function, q is a vector in R d , T d denotes the domain [-π, π] d with periodic boundary conditions, d = 1, 2 is the dimension, f (u) = u 2 /2, 0 < α ≤ 2 and (-∆) α/2 = Λ α is the fractional Laplacian.

This system was proposed by Othmers & Stevens [START_REF] Stevens | Aggregation, blowup, and collapse: the ABC's of taxis in reinforced random walks[END_REF] based on biological considerations as a model of tumor angiogenesis. In particular, in the previous system, u is the density of vascular endothelial cells and q = ∇ log(v) where v is the concentration of the signal protein known as vascular endothelial growth factor (VEGF) (see Bellomo, Li, & Maini [START_REF] Bellomo | On the foundations of cancer modelling: selected topics, speculations, and perspectives[END_REF] for more details on tumor modelling). Similar hyperbolic-dissipative systems arise also in the study of compressible viscous fluids or magnetohydrodynamics (see S. Kawashima [START_REF] Kawashima | Large-time behavior of solutions for hyperbolic-parabolic systems of conservation laws[END_REF] and the references therein).

Equation (1) appears as a singular limit of the following Keller-Segel model of aggregation of the slime mold Dictyostelium discoideum [START_REF] Keller | Initiation of slime mold aggregation viewed as an instability[END_REF] (see also Patlak [START_REF] Patlak | Random walk with persistence and external bias[END_REF]) [START_REF] Fleischer | Burgers' turbulence with self-consistently evolved pressure[END_REF] ∂ t u = -Λ α u -χ∇ • (u∇G(v)),

∂ t v = ν∆v + (f (u) + λ)v,
when G(v) = log(v) and the diffusion of the chemical is negligible, i.e. ν → 0. Similar equations arising in different context are the Majda-Biello model of Rossby waves [START_REF] Majda | The nonlinear interaction of barotropic and equatorial baroclinic Rossby waves[END_REF] or the magnetohydrodynamic-Burgers system proposed by Fleischer & Diamond [START_REF] Fleischer | Burgers' turbulence with self-consistently evolved pressure[END_REF].

Most of the results for (1) corresponds to the case where d = 1. Then, when the diffusion is local i.e. α = 2, (1) has been studied by many different research groups. In particular, Fan & Zhao [START_REF] Fan | Blow up criterion for a hyperbolic-parabolic system arising from chemotaxis[END_REF], Li & Zhao [START_REF] Li | Initial-boundary value problems for a system of hyperbolic balance laws arising from chemotaxis[END_REF], Mei, Peng & Wang [START_REF] Mei | Asymptotic profile of a parabolic-hyperbolic system with boundary effect arising from tumor angiogenesis[END_REF], Li, Pan & Zhao [START_REF] Li | Quantitative decay of a one-dimensional hybrid chemotaxis model with large data[END_REF], Jun, Jixiong, Huijiang & Changjiang [START_REF] Jun | Global solutions to a hyperbolic-parabolic coupled system with large initial data[END_REF] Li & Wang [START_REF] Li | Nonlinear stability of traveling waves to a hyperbolic-parabolic system modeling chemotaxis[END_REF] and Zhang & Zhu [START_REF] Zhang | Global existence of solutions to a hyperbolic-parabolic system[END_REF] studied the system (1) when α = 2 and f (u) = u under different boundary conditions (see also the works by Jin, Li & Wang [START_REF] Jin | Asymptotic stability of traveling waves of a chemotaxis model with singular sensitivity[END_REF], Li, Li & Wang [START_REF] Li | Stability of traveling waves of the Keller-Segel system with logarithmic sensitivity[END_REF], Wang & Hillen [START_REF] Wang | Shock formation in a chemotaxis model[END_REF] and Wang, Xiang & Yu [START_REF] Wang | Asymptotic dynamics on a singular chemotaxis system modeling onset of tumor angiogenesis[END_REF]). The case with general f (u) was studied by Zhang, Tan & Sun [START_REF] Zhang | Global existence and asymptotic behavior of smooth solutions to a coupled hyperbolic-parabolic system[END_REF] and Li & Wang [START_REF] Li | Nonlinear stability of large amplitude viscous shock waves of a generalized hyperbolic-parabolic system arising in chemotaxis[END_REF].

Equation [START_REF] Bellomo | On the foundations of cancer modelling: selected topics, speculations, and perspectives[END_REF] in several dimensions has been studied by Li, Li & Zhao [START_REF] Li | On a hyperbolic-parabolic system modeling chemotaxis[END_REF], Hau [START_REF] Hao | Global well-posedness for a multidimensional chemotaxis model in critical Besov spaces[END_REF] and Li, Pan & Zhao [START_REF] Li | Global dynamics of a hyperbolic-parabolic model arising from chemotaxis[END_REF]. There, among other results, the global existence for small initial data in H s , s > 2 is proved.

To the best of our knowledge, the only result when the diffusion is nonlocal, i.e. 0 < α < 2, is [START_REF] Granero-Belinchón | On the fractional fisher information with applications to a hyperbolic-parabolic system of chemotaxis[END_REF]. In that paper we obtained appropriate lower bounds for the fractional Fisher information and, among other results, we proved the global existence of weak solution for f (u) = u r /r and 1 ≤ r ≤ 2.

In this note, we address the existence of classical solutions in the case 0 < α ≤ 2. This is a challenging issue due to the hyperbolic character of the equation for q. In particular, u verifies a transport equation where the velocity q is one derivative more singular than u (so ∇•(uq) is two derivatives less regular than u).

Statement of the results

For the sake of clarity, let us state some notation: we define the mean as

g = 1 (2π) d T d g(x)dx.
Also, from this point onwards, we write H s for the L 2 -based Sobolev space of order s endowed with the norm

u 2 H s = u 2 L 2 + u 2 Ḣs , u Ḣs = Λ s u L 2 .
For β ≥ 0, we consider the following energies E β and dissipations D β , (4)

E β (t) = u 2 Ḣβ + q 2 Ḣβ , D β (t) = u 2 Ḣβ+α/2 .
Recall that the lower order norms verify the following energy balance [START_REF] Granero-Belinchón | On the fractional fisher information with applications to a hyperbolic-parabolic system of chemotaxis[END_REF] (5)

1 2 u(t) 2 L 2 + q(t) 2 L 2 + t 0 u(s) 2 Ḣα/2 ds = 1 2 u 0 2 L 2 + q 0 2 L 2 .
2.1. On the scaling invariance. Notice that the equations ( 1)-( 2) verify the following scaling symmetry: for every λ > 0

u λ (x, t) = λ α-1 u (λx, λ α t) , q λ (x, t) = λ α-1 q (λx, λ α t) .
This scaling serves as a zoom in towards the small scales. We also know that

u(t) 2 L 2 + q(t) 2 L 2
is the strongest (known) quantity verifying a global-in-time bound. Then, in the one dimensional case, the L 2 norms of u and q are invariant under the scaling of the equations when α = 1.5. That makes α = 1.5 the critical exponent for the global estimates known. Equivalently, if we define the rescaled (according to the scaling of the strongest conserved quantity

u(t) 2 L 2 + q(t) 2 L 2 ) functions u γ (x, t) = γ 0.5 u (γx, γ α t) , q γ (x, t) = γ 0.5 q (γx, γ α t) .
we have that u γ and q γ solve

∂ t u γ = -Λ α u γ + γ α-1.5 ∂ x (u γ q γ ), ∂ t q γ = γ α-1.5 u γ ∂ x u γ .
Larger values of α form the subcritical regime where the diffusion dominates the drift in small scales. Smaller values of α form the supercritical regime where the drift might be dominant at small scales.

Similarly, the two dimensional case has critical exponent α = 2.

Remark 1. Notice that the equations (1)-( 2) where f (u) = u have a different scaling symmetry but the same critical exponent α = 1.5. In this case, the scaling symmetry is given by

u λ (x, t) = λ 2α-2 u (λx, λ α t) , q λ (x, t) = λ α-1 q (λx, λ α t) ,
while the conserved quantity is u(t) L 1 + q(t) 2 L 2 /2. Thus, if we define the rescaled (according to the scaling of the conserved quantity) functions

u γ (x, t) = γu (γx, γ α t) , q γ (x, t) = γ 0.5 q (γx, γ α t) .
we have that u γ and q γ solve

∂ t u γ = -Λ α u γ + γ α-1.5 ∂ x (u γ q γ ), ∂ t q γ = γ α-1.5 ∂ x u γ .
A global existence result when α is the range 1.5 ≤ α < 2 for the problem where f (u) = u is left for future research.

Results in the one-dimensional case d = 1. One of our main results is

Theorem 1. Fix T an arbitrary parameter and let (u 0 , q 0 ) ∈ H 2 (T)×H 2 (T) be the initial data such that 0 ≤ u 0 and q 0 = 0. Assume that α ≥ 1.5. Then there exist a unique global solution (u(t), q(t)) to problem (1) verifying

u ∈ L ∞ (0, T ; H 2 (T)) ∩ L 2 (0, T ; H 2+α/2 (T)), q ∈ L ∞ (0, T ; H 2 (T)).
Furthermore, the solution is uniformly bounded in

(u, q) ∈ C([0, ∞), H 1 (T)) × C([0, ∞), H 1 (T)).
In the case where the strength of the diffusion, α, is even weaker, we have the following global existence result for small data: Theorem 2. Fix T an arbitrary parameter and let (u 0 , q 0 ) ∈ H 2 (T)×H 2 (T) be the initial data such that 0 ≤ u 0 and q 0 = 0. There exists C α such that if 1.5 > α > 1 and u 0 2 Ḣα/2 + q 0 2 Ḣα/2 ≤ C α then there exist a unique global solution (u(t), q(t)) to problem (1) verifying

u ∈ L ∞ (0, T ; H 2 (T)) ∩ L 2 (0, T ; H 2+α/2 (T)), q ∈ L ∞ (0, T ; H 2 (T)).
Furthermore, the solution verifies

u(t) 2 Ḣα/2 + q(t) 2 Ḣα/2 ≤ u 0 2 Ḣα/2 + q 0 2 Ḣα/2 .
Corollary 1. Fix T an arbitrary parameter and let (u 0 , q 0 ) ∈ H 2 (T)×H 2 (T) be the initial data such that 0 ≤ u 0 and q 0 = 0. Assume that 1 ≥ α ≥ 0.5 and

u 0 2 Ḣ1 + q 0 2 Ḣ1 < 4 9C 2 S
where C S is defined in [START_REF] Jun | Global solutions to a hyperbolic-parabolic coupled system with large initial data[END_REF]. Then there exist a unique global solution (u(t), q(t)) to problem (1) verifying

u ∈ L ∞ (0, T ; H 2 (T)) ∩ L 2 (0, T ; H 2+α/2 (T)), q ∈ L ∞ (0, T ; H 2 (T)).
Furthermore, the solution verifies

u(t) 2 H 1 + q(t) 2 H 1 ≤ u 0 2 H 1 + q 0 2 H 1 . 2.3.
Results in the two-dimensional case d = 2. In two dimensions the global existence read Theorem 3. Fix T an arbitrary parameter and let (u 0 , q 0 ) ∈ H 2 (T 2 ) × H 2 (T 2 ) be the initial data such that 0 ≤ u 0 , q 0 = 0 and curl q 0 = 0. Assume that α = 2. Then there exist a unique global solution (u(t), q(t)) to problem (1) verifying

u ∈ L ∞ (0, T ; H 2 (T 2 )) ∩ L 2 (0, T ; H 3 (T 2 )), q ∈ L ∞ (0, T ; H 2 (T 2 )).
Furthermore, the solution is uniformly bounded in

(u, q) ∈ C([0, ∞), H 1 (T 2 )) × C([0, ∞), H 1 (T 2 )).
Corollary 2. Fix T an arbitrary parameter and let (u 0 , q 0 ) ∈ H 2 (T 2 ) × H 2 (T 2 ) be the initial data such that 0 ≤ u 0 , q 0 = 0 and curl q 0 = 0. Assume that 2 > α ≥ 1 and

u 0 2 Ḣ1 + q 0 2 Ḣ1 < C
where C is a universal constant. Then there exist a unique global solution (u(t), q(t)) to problem (1) verifying

u ∈ L ∞ (0, T ; H 2 (T 2 )) ∩ L 2 (0, T ; H 2+α/2 (T 2 )), q ∈ L ∞ (0, T ; H 2 (T 2 )).
Furthermore, the solution verifies

u(t) 2 H 1 + ∇ • q(t) 2 L 2 ≤ u 0 2 H 1 + ∇ • q 0 2 L 2 .
Remark 2. In the case where the domain is the one-dimensional torus, T, local existence of solution for (1)-( 2) was proved in [START_REF] Granero-Belinchón | On the fractional fisher information with applications to a hyperbolic-parabolic system of chemotaxis[END_REF] for a more general class of kinetic function f (u). The local existence of solution for (1)-( 2) the domain is the two-dimensional torus T d with d = 2 follows from the local existence result in [START_REF] Granero-Belinchón | On the fractional fisher information with applications to a hyperbolic-parabolic system of chemotaxis[END_REF] with minor modifications. Consequently, we will focus on obtaining global-in-time a priori estimates.

2.4. Discussion. Due to the hyperbolic character of the equation for q, prior available global existence results of classical solution for equation ( 1) impose several assumptions. Namely,

• either d = 1 and α = 2 [START_REF] Zhang | Global existence and asymptotic behavior of smooth solutions to a coupled hyperbolic-parabolic system[END_REF][START_REF] Li | Nonlinear stability of large amplitude viscous shock waves of a generalized hyperbolic-parabolic system arising in chemotaxis[END_REF],

• or d = 2, 3, α = 2 and the initial data verifies some smallness condition on strong Sobolev spaces H s , s ≥ 2 [START_REF] Xie | Global existence and convergence rates for the strong solutions into the 3d chemotaxis model[END_REF][START_REF] Zhang | Global existence and exponential stability for the strong solutions in H 2 to the 3-d chemotaxis model[END_REF]. Our results removed some of the previous conditions. On the one hand, we prove global existence for arbitrary data in the cases d = 1 and α ≥ 1.5 and d = 2 and α = 2. On the other hand, in the cases where we have to impose size restrictions on the initial data, the Sobolev spaces are bigger than H 2 (thus, the norm is weaker). Finally, let us emphasize that our results can be adapted to the case where the spatial domain is R d .

A question that remains open is the trend to equilibrium. From ( 5) is clear that the solution (u(t), q(t)) tends to the homogeneous state, namely ( u 0 , 0). However, the rate of this convergence is not clear.

Proof of Theorem 1

Step 1; H 1 estimate: Testing the first equation in (1) against Λ 2 u, integrating by parts and using the equation for q, we obtain 1 2

d dt u 2 Ḣ1 + ∂ x u 2 Ḣα/2 = - T ∂ x (uq)∂ 2 x udx = 1 2 T ∂ x q(∂ x u) 2 dx - T ∂ x qu∂ 2 x udx = 1 2 T ∂ x q(∂ x u) 2 dx - T ∂ x q(∂ t ∂ x q -(∂ x u) 2 )dx, so 1 2 d dt ( u 2 Ḣ1 + q 2 Ḣ1 ) + u 2 Ḣ1+α/2 = 3 2 T ∂ x q(∂ x u) 2 dx. Denoting I = 3 2 T ∂ x q(∂ x u) 2 dx,
and using Sobolev embedding and interpolation, we have that (6)

I ≤ 3 2 q Ḣ1 ∂ x u 2 L 4 ≤ 3 2 C S q Ḣ1 ∂ x u 2 Ḣ0.25 ,
where C S is the constant appearing in the embedding [START_REF] Jun | Global solutions to a hyperbolic-parabolic coupled system with large initial data[END_REF] g L 4 ≤ C S g Ḣ0. [START_REF] Zhang | Global existence of solutions to a hyperbolic-parabolic system[END_REF] .

Using the interpolation

H 1+α/2 ⊂ H 1.25 ⊂ H α/2 ,
and Poincaré inequality (if α > 1.5) we conclude

I ≤ c q Ḣ1 Λ α/2 u L 2 u Ḣ1+α/2 ,
Using (4), we have that

d dt E 1 + D 1 ≤ c u 2 Ḣα/2 E 1 .
Using Gronwall's inequality and the estimate ( 5), we have that

sup 0≤t<∞ E 1 (t) ≤ C( u 0 H 1 , q 0 H 1 ), T 0 D 1 (s)ds ≤ C( u 0 H 1 , q 0 H 1 , T ), ∀0 < T < ∞.
Step 2; H 2 estimate: Now we prove that the solutions satisfying the previous bounds for E 1 and D 1 also satisfy the corresponding estimate in H 2 . We test the equation for u against Λ 4 u. We have that 1 2

d dt u 2 Ḣ2 + u 2 Ḣ2+α/2 = - T ∂ 2 x (uq)∂ 3 x udx = T ∂ x q 5(∂ 2 x u) 2 2 dx - T ∂ 2 x q(∂ t ∂ 2 x q -5∂ x u∂ 2 x u)dx, so 1 2 d dt ( u 2 Ḣ2 + q 2 Ḣ2 ) + u 2 Ḣ2+α/2 = 5 2 T ∂ x q(∂ 2 x u) 2 dx + 5 T ∂ 2 x q∂ 2 x u∂ x udx.
We define

J 1 = 5 2 T ∂ x q(∂ 2 x u) 2 dx, J 2 = 5 T ∂ 2 x q∂ 2 x u∂ x udx.
Then, we have that

J 1 ≤ c ∂ x q L ∞ ∂ 2 x u 2 L 2 ≤ c ∂ 2 x q 0.5 L 2 u α Ḣ1+α/2 u 2-α Ḣ2+α/2 ,
so, using Young's inequality,

J 1 ≤ c ∂ 2 x q 1 α L 2 u 2 Ḣ1+α/2 + 1 4 u 2 Ḣ2+α/2 .
Similarly, using Poincaré inequality and α ≥ 0.5,

J 2 ≤ c ∂ x u L 4 ∂ 2 x u L 4 ∂ 2 x q L 2 ≤ c u Ḣ1+α/2 u Ḣ2+α/2 ∂ 2 x q L 2 ,
and

J 2 ≤ c u 2 Ḣ1+α/2 ∂ 2 x q 2 L 2 + 1 4 u 2 Ḣ2+α/2 .
Finally,

d dt E 2 (t) + D 2 (t) ≤ c u 2 Ḣ1+α/2 (E 2 (t) + 1)
and we conclude using Gronwall's inequality.

Proof of Theorem 2

Step 1; H α/2 estimate: Testing the first equation in (1) against Λ α u, we obtain 1 2

d dt u 2 Ḣα/2 + u 2 Ḣα = T ∂ x (uq)Λ α udx = - T Λ α (uq)∂ x udx = - T (Λ α (uq) -uΛ α q) ∂ x udx - T Λ α qu∂ x udx, so 1 2 d dt u 2 Ḣα/2 + q 2 Ḣα/2 + u 2 Ḣα ≤ - T [Λ α , u]q∂ x udx.
We define

K = - T [Λ α , u]q∂ x udx.
Using the classical Kenig-Ponce-Vega commutator estimate [START_REF] Kenig | Well-posedness and scattering results for the generalized Korteweg-De Vries equation via the contraction principle[END_REF] and Sobolev embedding, we have that

[Λ α , u]q L 2 ≤ c ∂ x u L 2+ǫ Λ α-1 q L 4+2ǫ ǫ + Λ α u L 2 q L ∞ ≤ c u Ḣ1+ ǫ 4+2ǫ Λ α-1 q Ḣ 1 2 -ǫ 4+2ǫ + u Ḣα q Ḣ α 2 . (8)
Thus, taking ǫ such that 1 + ǫ 4 + 2ǫ = α, i.e. ǫ = 4α -4 3 -2α Equation ( 8) reads [START_REF] Keller | Initiation of slime mold aggregation viewed as an instability[END_REF] [Λ α , u]q L 2 ≤ c u Ḣα q Ḣ 1 2 + u Ḣα q Ḣ α 2 Using (9) and Poincaré inequality, we have that

K ≤ c u 2 Ḣα q Ḣ α 2
Then, we have that

d dt E α 2 + D α 2 ≤ c E α 2 D α 2 .
Thus, due to the smallness restriction on the initial data, we obtain

E α 2 (t) + δ t 0 D α 2 (s)ds ≤ E α 2 (
0) for 0 < δ small enough.

Step 2; H 1 estimate: Our starting point is [START_REF] Jin | Asymptotic stability of traveling waves of a chemotaxis model with singular sensitivity[END_REF]. Then we use the interpolation g 2 Ḣ0.25 ≤ c g L 2 g Ḣ0. [START_REF] Hao | Global well-posedness for a multidimensional chemotaxis model in critical Besov spaces[END_REF] , to obtain

I ≤ c q Ḣ1 u Ḣ1 u Ḣ1.5 ≤ cE 1 D α 2 + D 1 2 .
Collecting all the estimates, we have that

d dt E 1 + D 1 ≤ cE 1 D α 2 ,
and we conclude using Gronwall's inequality. The H 2 estimates follows as in the proof of Theorem 1.

Proof of Corollary 1

Using α ≥ 0.5 and the estimate (6), we have that

I ≤ 3 2 C S q Ḣ1 u 2 Ḣ1+α/2 ≤ 3 2 C S E 1 D 1 .
Thus, 1 2

d dt E 1 + D 1 ≤ 3 2 C S E 1 D 1 .
Thus, due to the smallness restriction on the initial data, we obtain

E 1 (t) + δ t 0 D 1 (s)ds ≤ E 1 (0)
for 0 < δ small enough. Equipped with this estimates, we can repeat the argument as in Step 2 in Theorem 1.

Proof of Theorem 3

Recall that the condition curl q 0 = 0 propagates in time, i.e. curl q(t) = curl q 0 + 1 2 t 0 curl∇u 2 ds = 0.

Using Plancherel Theorem, we have that

∇q 2 L 2 = C ξ∈Z 2 |ξ| 2 |q(ξ)| 2 = C ξ∈Z 2 (ξ 2 1 + ξ 2 2 )(q 2 1 + q2 2 ).
Due to the irrotationality ξ ⊥ • q = 0. Then, we compute

∇ • q 2 L 2 = C ξ∈Z 2 |ξ • q(ξ)| 2 = C ξ∈Z 2 (ξ 1 q1 (ξ) + ξ 2 q2 (ξ)) 2 = C ξ∈Z 2 (ξ 1 q1 (ξ)) 2 + (ξ 2 q2 (ξ)) 2 + 2ξ 1 q1 (ξ)ξ 2 q2 (ξ) = C ξ∈Z 2 (ξ 1 q1 (ξ)) 2 + (ξ 2 q2 (ξ)) 2 + (ξ 2 q1 (ξ)) 2 + (ξ 1 q2 (ξ)) 2 .
So, the vector field q satisfies

∇q L 2 ≤ ∇ • q L 2 .
As a consequence of ∂ t q i = 0 and q 0 = 0, every coordinate of q satisfy q i (t) = 0, and the Poincaré-type inequality

(10) q L 2 ≤ c ∇ • q L 2 .
Notice that in two dimensions we also have the energy balance [START_REF] Hao | Global well-posedness for a multidimensional chemotaxis model in critical Besov spaces[END_REF]. We test equation (1) against Λ 2 u and use the equation for q. We obtain 1 2

d dt u 2 Ḣ1 + ∇ • q 2 L 2 = -u 2 Ḣ2 - T 2 ∇u • q∆udx + T 2 |∇u| 2 ∇ • qdx.
Using Hölder inequality, Sobolev embedding and interpolation, we have that

d dt u 2 Ḣ1 + ∇ • q 2 L 2 + 2 u 2 Ḣ2 ≤ c u Ḣ1.5 q L 4 u Ḣ2 + u 2 Ḣ1.5 ∇ • q L 2 ≤ c u 0.5 Ḣ1 q 0.5 L 2 q 0.5 H 1 u 1.5 Ḣ2 + c u Ḣ1 u Ḣ2 ∇ • q L 2 .
Using the Hödge decomposition estimate together with the irrotationality of q and (10), we have that ( 11)

q H 1 ≤ c ( q L 2 + ∇ • q L 2 ) ≤ c ∇ • q L 2 .
Due to (5), we obtain that

q 2 L ∞ (0,∞,L 2 ) + u 2 L 2 (0,∞, Ḣ1 ) ≤ C( u 0 L 2 , q 0 L 2 ) so, d dt u 2 Ḣ1 + ∇ • q 2 L 2 + u 2 Ḣ2 ≤ c u 2 Ḣ1 ∇ • q 2 L 2 .
Using Gronwall's inequality and the integrability of u 2 Ḣ1 (see ( 5)), we obtain E 1 ≤ C( u 0 H 1 , q 0 H 1 ), T 0 D 1 (s)ds ≤ C( u 0 H 1 , q 0 H 1 , T ), ∀0 < T < ∞.

To obtain the H 2 estimates, we test against Λ 4 u. Then, using the previous H 1 uniform bound and q 2 L ∞ ≤ c q L 2 q H 2 ≤ c q L 2 ∆q L 2 , we have that [START_REF] Li | Quantitative decay of a one-dimensional hybrid chemotaxis model with large data[END_REF] 1 2

d dt ∆u 2 L 2 + u 2 Ḣ3 = - T d ∇∆u∇(∇u • q)dx - T d u∇∆u • ∇(∇ • q)dx - T d ∇u • ∇∆u∇ • qdx.
Due to the irrotationality of q and the identity ∇∇ • q -∆q = curl (curl q), we have ∂ t ∆q = ∂ t ∇(∇ • q) = ∇|∇u| 2 + u∇∆u + ∇u∆u. Applying Sobolev embedding and interpolation, we obtain that [START_REF] Li | Quantitative decay of a one-dimensional hybrid chemotaxis model with large data[END_REF] and we conclude using Gronwall's inequality.

Proof of Corollary 2

We test the equation (1) against Λ 2 u. We obtain that 1 2

d dt u 2 Ḣ1 + ∇ • q 2 L 2 = -u 2 Ḣ1+ α 2 + T 2 ∇(∇u•q)∇udx+ T 2 |∇u| 2 ∇•qdx.
After a short computation, using Hölder estimates, Sobolev embedding and interpolation, we obtain that 1 2

d dt u 2 Ḣ1 + ∇ • q 2 L 2 + u 2 Ḣ1+ α 2 ≤ c u 2 Ḣ1.5 q Ḣ1 .
Using [START_REF] Li | On a hyperbolic-parabolic system modeling chemotaxis[END_REF] and α ≥ 0, we obtain 1 2

d dt u 2 Ḣ1 + ∇ • q 2 L 2 + D 1 ≤ cD 1 ∇ • q 2 L 2 .
We conclude the result with the previous ideas.

  + D 2 ≤ cE 2 ,

	so,
		d dt	E 2
			can be
	estimated as
	d dt	∆u 2 L 2 + ∆q 2 L 2 + u 2 Ḣ3 ≤ c ∆q 2 L 2 ,

Acknowledgment

The author is funded by the Labex MILYON and the Grant MTM2014-59488-P from the Ministerio de Economía y Competitividad (MINECO, Spain).


Leandro R Lustosa
email: leandro-ribeiro.lustosa@isae.fr

Sergio Pizziol
email: sergio.pizziol@isae.fr

François Defaÿ
email: francois.defay@isae.fr

Jean-Marc Moschetta
email: jean-marc.moschetta@isae.fr

An Error Model of a Complementary Filter for use in Bayesian Estimation -The CF-EKF Filter

Keywords: Sensor data fusion, filtering, estimation, UAV navigation, guidance and control

It is well known that stand-alone inertial navigation systems (INS) have their errors diverging with time. Consequently, an upper bound on the duration of INS systems precludes their use in low-cost micro unmanned aerial vehicles. The traditional approach for solving this matter is to resort to aiding devices, e.g., global navigation satellite system (GNSS) receivers, sighting devices, etc. Two philosophies have been extensively applied to perform data fusion: extended Kalman filtering (EKF) and complementary filtering (CF). Previous work in the literature showed that the computationally less expensive CF can be robustly applied to attitude estimation using low-cost sensors and achieve performance that is comparable to that of a full EKF. However, performance is degraded by vehicle manoeuvres and no measurement on estimate uncertainties is given. Furthermore, a large number of sensors makes it impracticable for optimal tuning of the CF. The present work lays the foundation for sensor filtering that employs the CF for attitude estimation by means of a magnetometer as an external aid, and an EKF for additional sensors integration. The main feature of this architecture is the possibility of deployment in a distributed multi-platform system and implementation of fault isolation by running the CF stage in a separate low-throughput reliable machine for stand-alone degraded mode operation. A case study is performed on synthetic data from inertial, magnetic and GNSS sensors.

INTRODUCTION

The dead-reckoning nature of stand-alone inertial navigation systems (INS) impose cumulative errors in the estimation of position, velocity and attitude. Consequently, an upper bound on the duration of stand-alone INS systems precludes their use in low-cost or micro unmanned aerial vehicles (µUAV). The traditional approach for solving such inconvenience is to resort to aiding devices such as global navigation satellite system (GNSS) receivers, magnetometers and/or sighting devices. Two philosophies have been extensively applied to perform this data fusion: extended Kalman filtering (EKF) and complementary filtering (CF).

Some key similarities and differences between EKF and CF are discussed in [START_REF] Higgins | A comparison of complementary and Kalman filtering[END_REF]]. For instance, a CF is equivalent to an EKF under steady-state conditions. However, navigation errors dynamics and observability are dependent on vehicle trajectory [START_REF] Goshen-Meskin | Observability analysis of piece-wise constant systems, ii. application to inertial navigation in-flight alignment[END_REF]] in such a way that an EKF can outperform a CF in manoeuvring conditions. Be that as it may, the use of low-cost sensors renders the EKF suboptimal and often inconsistent due to large navigation error covariances that undermine linearization assumptions. Particle filters [START_REF] Gordon | Novel approach to nonlinear/non-Gaussian Bayesian state estimation[END_REF]] have been employed to solve such This work is partially supported by the Conselho Nacional de Desenvolvimento Científico e Tecnológico, CNPq.

problem but require a fair amount of computational power often not available in low-cost or micro-UAVs applications. CF, on the other hand, is an exceptionally inexpensive filter.

Previous work [START_REF] Euston | A complementary filter for attitude estimation of a fixed-wing uav[END_REF]] has shown that the CF can be robustly applied to low-cost sensors and achieve performance that is comparable to that of a full suboptimal EKF based on GNSS-aided inertial measurement unit (IMU/GNSS) when vehicle dynamics is taken into consideration by means of airspeed measurements. Another strategy to relax the CF required assumption of constant velocity motion is to apply gains that adapt to the motion to be estimated [START_REF] Calusdian | Adaptive-gain complementary filter of inertial and magnetic data for orientation estimation[END_REF]]. The issue of tuning a given CF is discussed in detail in [START_REF] Vasconcelos | Discrete time-varying attitude complementary filter[END_REF]] to shape a frequency response that blends the frequency contents of the aiding devices and the inertial sensors (implementation aspects are also highlighted). [START_REF] Mahony | Complementary filter design on the special orthogonal group so(3)[END_REF]] compares two implementation architectures for the CF, namely, direct CF and passive CF, and shows that the latter has superior dynamical properties, proves its convergence and extends it to provide adaptive rate-gyro drift estimation.

The present work exploits the passive CF for position, velocity and attitude estimation employing a magnetometer as an external aid to bound attitude errors. On top of that, an external EKF is applied to restrain position and velocity errors by means of GNSS measurements. To the authors' knowledge, there is no previous work on how to seamlessly connect both filters, although there are ad hoc CF-EKF filter interconnection and tuning [START_REF] Jung | Inertial attitude and position reference system development for a small uav[END_REF]]. This paper contributes by delineating the overall CF-EKF filter structure designed to simplify the tedious tuning process, studying some of its properties and delivering implementation details. Furthermore, some EKF aiding devices (sighting devices, for instance) call for operational systems that increases system complexity thus reducing INS reliability. In view of this, it is shown how to effectively employ CF-EKF in a multi-platform UAV to implement fault isolation by running the CF stage in a separate low-throughput low-level reliable machine for stand-alone degraded mode operation.

In the following, an attitude passive CF based on gravitational and magnetic fields is reviewed (section 3.1) and its dynamic properties are discussed (section 3.3). Moreover, this work contributes by delivering an error model of this complementary filter based on psi-angle representation (section 3.2) for statistical error analysis and extended Kalman or information filter integration (section 4). The proposed framework is evaluated by means of a distributed loosely-coupled IMU/GNSS Kalman filter applied to synthetic IMU/GNSS sensor data (section 5). Finally, concluding remarks and perspectives are presented (section 6).

NOTATION AND REFERENCE FRAMES

In this paper, kinematic quantities of interest in multiple moving reference frames are algebraically studied and, accordingly, a consistent and precise notation is called for. Therefore, this section states the notation conventions applied herein, including relevant reference frames definition and related variables of interest.

The notation a x c is employed, where the symbol x is replaced by the symbol of the desired vector quantity (p for position, v for velocity, a for acceleration, ω for angular velocity) of frame/point (depending on the context) C with respect to frame A. For instance, i ω b denotes angular velocity of frame B with respect to frame I.

The decomposition of a vector x ∈ R n into its components in a coordinate system R is denoted by means of the right subscript position, e.g.

x r = (x r1 x r2 • • • x rn ) T (1)
We make extensive use of the vector product operation and its matrix representation (in some basis B) is denoted by

v b × = [v b ×] = 0 -v b3 v b2 v b3 0 -v b1 -v b2 v b1 0 (2)
Quaternion algebra structure definition varies to a small extent in literature, but herein a quaternion q ∈ (R 4 , ×) is defined as [START_REF] Stevens | Aircraft control and simulation[END_REF]]

q = q 0 q (3)
where q 0 ∈ R and q ∈ R 3 with quaternion product operation defined as p × q = p 0 q 0 -p • q p 0 q + q 0 p + p × q (4)

The quaternion and direction cosine matrix (DCM) representations of rotation from frame A to B, respectively q a b and D a b , are defined in the direction such that the rotation formulas [START_REF] Stevens | Aircraft control and simulation[END_REF]] are written as The inertial frame I is assumed stationary with respect to Earth since this work targets low-cost sensors that possess low rate-gyro precision when compared to Earth angular velocity order-of-magnitude. On the other hand, the body axes B are fixed in the inertial measurement unit (IMU) and yields a natural vector space basis for sensor measurements. Finally, an embedded processor collects sensor data and computes IMU position and attitude with respect to I by means of the attitude complementary filter (see section 3). Due to sensors errors, the estimated orientation of B is misaligned and it is denoted by the computed frame C.

0 x b = (q a b ) × 0 x a × q a b (

ATTITUDE COMPLEMENTARY FILTER

Filter description

This section reviews an attitude complementary filter that, when applied to rate-gyro, accelerometer and magnetometer measurements, yields attitude estimation with bounded errors (see figure 2). Additionally, it delivers position and velocity estimators. However, their errors grow with time and call for additional sensors and data fusion algorithms (see section 4).

Complementary

Filter This work models the rate-gyro measurement ω(t) as a superposition of the nominal value i ω b b (t), a random constant (with respect to body axes) Gaussian drift ε b ∼ N (0, Σ ε ) and white Gaussian noise

ω(t) f (t) b(t) q i c (t) p c i (t) i v c i (t)
ν g b (t) ∼ N (0, Σ g ) such that ω(t) = i ω b b (t) + ε b + ν g b (t) (7)
Similarly, accelerometer and magnetometer measurements, f (t) and b(t), respectively, are modeled as

f (t) = i a b b (t) -g b (t) + ∇ b + ν a b (t) (8) and b(t) = B b (t) + ∆ b + ν m b (t) (9) where ∇ b ∼ N (0, Σ ∇ ), ν a b (t) ∼ N (0, Σ a ), ∆ b ∼ N (0, Σ ∆ ), ν m b (t) ∼ N (0, Σ m )
, and g denote, respectively, accelerometer/magnetometer bias/noise, and Earth gravitational acceleration.

Notice that direct integration of equation 7 yields attitude estimators with unbounded errors that call for additional sensors. The present version of the attitude complementary filter assumes the local magnetic field B i is known and compares a predicted measure B c with an ideal measure B b to produce a residual that should be closely related to attitude errors. Mathematically, notice that

B b × B c = -[B c ×]B b = -[B c ×]D c b B c (10) where D c
b is the DCM that transforms computed frame coordinates into body axes components. Assuming a 3 -2 -1 rotation order of ψ 3 , ψ 2 and ψ 1 , respectively, D c b can be written as

D c b = cψ 2 cψ 3 cψ 2 sψ 3 -sψ 2 sψ 1 sψ 2 cψ 3 -cψ 1 sψ 3 cψ 1 cψ 3 + sψ 1 sψ 2 sψ 3 sψ 1 cψ 2 sψ 1 sψ 3 + cψ 1 sψ 2 cψ 3 cψ 1 sψ 2 sψ 3 -sψ 1 cψ 3 cψ 1 cψ 2 (11)
If small error angles ψ i are assumed, cosψ i → 1, sinψ i → ψ i , and equation 11 is simplified to

D c b = 1 ψ 3 -ψ 2 -ψ 3 1 ψ 1 ψ 2 -ψ 1 1 = I -[ψ×] ( 12
)
where

ψ = (ψ 1 ψ 2 ψ 3 ) T (13)
Therefore, equation 10 can be simplified to read

B b × B c = -[B c ×](I -[ψ×])B c = -[B c ×] 2 ψ (14)
Assuming the local gravitational vector is known, a similar development yields

g b × g c = -[g c ×] 2 ψ (15)
Equations 14 and 15 contain candidates for feedback stabilization of ψ as a consequence of their proportionality to ψ. Accordingly, an estimator i ω c c for i ω b b , namely, the angular velocity of the body reference frame with respect to the inertial frame described in the body axes, is defined as

i ω c c = ω -k a f × g c + k m b × B c (16)
where k a , k m > 0. From the embedded processor point of view, equation 16 can be computed as

i ω c c = ω -k a f × D i c g i + k m b × D i c B i (17) where D i
c is computed by means of q i c , which can be numerically computed from the following classic [START_REF] Stevens | Aircraft control and simulation[END_REF]] quaternion differential equation

d dt q i c = 1 2 0 -( i ω c c ) T i ω c c -[ i ω c c ×] q i c ( 18
)
and quaternion to DCM transformation equation

D i c =   (1 -2q 2 2 -2q 2
3 ) 2(q 1 q 2 + q 0 q 3 ) 2(q 1 q 3 -q 0 q 2 ) 2(q 1 q 2 -q 0 q 3 ) (1 -2q 2 1 -2q 2 3 ) 2(q 2 q 3 + q 0 q 1 ) 2(q 1 q 3 + q 0 q 2 ) 2(q 2 q 3 -q 0 q 1 ) (1

-2q 2 1 -2q 2 2 )   (19)
Finally, position and velocity estimators are obtained by means of direct integration of the following mechanization equations d dt

p c i = i v c i ( 20
)
and

d dt i v c i = D c i f + g i (21)
Equations 17, 18, 19, 20 and 21 implement the CF inertial navigation algorithm of figure 2. In the next two subsections, some of its properties are re-derived in view of the misalignment model, which is the main contribution of this work and the necessary foundation for CF-EKF interfacing.

Complementary filter misalignment error model

This section develops an error model for ψ(t), i.e., a stochastic differential equation that models complementary filter inertial navigation errors in time due to vehicle motion and uncertainties in sensors.

Initially, notice that ψ(t) is composed by Euler angles and their relation to angular velocity is dictated by [START_REF] Stevens | Aircraft control and simulation[END_REF]]

d dt ψ =    1 tan ψ 2 sin ψ 1 tan ψ 2 cos ψ 1 0 cos ψ 1 -sin ψ 1 0 sin ψ 1 cos ψ 2 cos ψ 1 cos ψ 2    c ω b b (22)
Assuming small angles and gathering only first order terms, we can write

d dt ψ = c ω b b (23) so that d dt ψ = i ω b b -i ω c b = i ω b b -D c b • i ω c c = i ω b b -i ω c c + [ψ×] i ω c c
(24) by means of equation 12. Therefore,

i ω b b -i ω c c = d dt ψ + [ i ω c c ×]ψ (25)
Substituting equations 7, 8, 16 into 25, we obtain

dψ dt = -[ i ω c c ×] + k m [B c ×] 2 + k a [g c ×] 2 ψ+ + -I ε b + -k a [g c ×] ∇ b + k m [B c ×] ∆ b + + -I ν g b + k m [B c ×] ν m b + -k a [g c ×] ν a b + + -k a [g c ×] i a b b (26)
On the other hand, the errors on position and velocity, namely, δp p cp b (27) and δv i v ci v b (28) have the following well-known differential equations for their time evolution [START_REF] Bar-Shalom | Estimation with aapplication to tracking and navigation: theory, algorithms and software[END_REF]]

d dt δp i = δv i (29) and d dt δv i = D c i [f ×]ψ + D c i ∇ b + D c i ν a b ( 30
)
Equations 26, 29 and 30 describe errors stochastic time evolution in the CF inertial navigation algorithm and can be used for several purposes, e.g., algorithm stability analysis (subsection 3.3) and Kalman filtering data fusion (section 4).

Convergence analysis

Given sensors noise and bias characteristics, equation 26 provides the means to statistically predict the performance of the CF (for small ψ). Furthermore, it describes the dynamics of the attitude estimation error ψ(t). Considering bias and noise quantities as filter inputs, the local stability of the attitude CF is determined by the following linear time-variant system

dψ dt = -[ i ω c c ×] + k m [B c ×] 2 + k a [g c ×] 2 N ψ (31)
Although ψ has no physical meaningful basis (among other, rotations do not have a vector space structure), one can assign a basis to ψ and treat it as if it was a vector algebraically. Since equation 31 is fully described in computed basis, we assign the C basis to ψ such that

dψ c dt + [ i ω c c ×]ψ c = k m [B c ×] 2 + k a [g c ×] 2 ψ c ( 32
)
The reasoning that transformed equation 31 into 32 is subtle but allows one to evoke the transport theorem from classical mechanics and rewrite the left-hand side as

dψ i dt c = k m [B c ×] 2 + k a [g c ×] 2 ψ c (33) such that dψ i dt = D c i k m [B c ×] 2 + k a [g c ×] 2 M D i c ψ i (34)
The interested reader can quickly check that every matrix of the form [v×] 2 has eigenvalues λ 1 = 0, with algebraic multiplicity 1, and λ 2 = -|v| 2 with algebraic multiplicity 2. Therefore, M is a linear combination of such matrices and it is negative semi-definite if k m > 0 and k a > 0. Therefore, the CF filter is marginally stable1 for attitude determination.

Position and velocity estimation errors, on the other hand, remains unbounded and call for additional sensors, which are commonly integrated by means of Kalman filters. The next section provides a clean strategy to easily integrate and tune an EKF on top of a complementary filter.

KALMAN FILTER INTEGRATION

In the light of the foregoing development, it is convenient to define the EKF state vector as (36) such that ẋ = Ax + Bw (37) where A and B can be obtained by inspection of equations 26, 29 and 30, yielding

x = (δp i δv i ψ ∇ b ε b ∆ b ) T (
A =        0 I 0 0 0 0 0 0 D c i [f ×] D c i 0 0 0 0 N -k a [g c ×] -I k m [B c ×] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0        (38)
and

B =        0 0 0 0 D c i 0 0 0 -k a [g c ×] -I k m [B c ×] -k a [g c ×] 0 0 0 0 0 0 0 0 0 0 0 0        (39)
On the other hand, the EKF observation equation is dependent on the available aiding sensors. In this paper, for illustration, a loosely-coupled GNSS receiver architecture [START_REF] Farrell | The global positioning system and inertial navigation[END_REF]] is employed such that y = Hx + v (40) where

H = I 3×3 0 3×3 0 3×12 0 3×3 I 3×3 0 3×12 (41)
and v is dependent on GNSS receiver uncertainties. After discretization of equation 37, the EKF is applied (for a comprehensive discussion on Kalman filtering and useful equivalent alternative formulations for use in low-cost architectures, see [START_REF] Maybeck | Stochastic models, estimation, and control[END_REF]]) and estimation of navigation errors and sensor biases is performed. After each EKF update step, the state x is used for CF correction (position, velocity and attitude feedback) and sensor on-line calibration. Figure 3 illustrates the overall architecture.

A remarkable feature of the CF-EKF is that given the statistics of the noise process defined by equation Notice that errors and covariances for position estimation are similar throughout all three implementations. In particular, the position covariance for the CF-EKF is only slightly larger than the optimal EKF. This feature illustrates the trade-off of centralized optimal and distributed sub-optimal estimation architectures. This tradeoff is stronger in velocity estimation where a factor of 2 is observed between EKF and CF-EKF estimated error covariances. Kalman filter statistical divergence (notice that RMS divergence is also observed). This precludes filter validation and usage. Suboptimal inflation of process noise could bring the filter back to statistical convergence, but the tuning strategy is often a blind and painful trial-and-error method which often yield poor overall performance. Similar performance patterns are observed in misalignment and accelerometer error curves.

Statistical inconsistency and RMS error divergence is also seen in ad hoc CF/EKF estimation of rate-gyro drift. Moreover, the trade-off in this component is observability: while this vehicle trajectory allows for drift observability in EKF, the CF-EKF does not react to it. Finally, none of the filters has magnetic bias observability.

CONCLUSION

This paper laid the foundation for CF and EKF interfacing without imposing excessive suboptimal assumptions that often reduce performance and make it impracticable for tuning. While CF employed IMU/magnetometer measurements to bound attitude errors, an external EKF bounded position and velocity errors by means of CF output and GNSS measurements. Any other traditional EKF aiding sensor could be employed in the CF-EKF architecture as well.

The main feature of the CF-EKF algorithm is the possibility of deployment in a distributed multi-platform system. This paper showed how to implement fault isolation by running the CF stage in a separate low-throughput low-level reliable machine for CF stand-alone degraded mode operation. This allows for manual recovery of an auto-pilot failure due to a navigation system fault. Although the CF is simple enough to be directly implemented in a simple computing hardware, the EKF calls for high-order matrices inversion that are numerically expensive and unstable. Ongoing work investigates equivalent numerically more robust approaches to the classic EKF formulation (such as U-D factorization) for implementation in an embedded µUAV real-time computer [START_REF] Lustosa | Longitudinal study of a tilt-body vehicle: modeling, control and stability analysis[END_REF]]. This will allow for proper assessment of its performance in more realistic conditions under real-time control system requirements (e.g., continuous position and velocity estimates).

Fig. 1 .

 1 Fig. 1. Reference frames illustration. The inertial (I), body (B) and computed (C) frames are fixed to the Earth, IMU and estimated orientation by the attitude complementary filter, respectively.

Fig. 2 .

 2 Fig. 2. Input-output schematic view of the quaternion complementary filter.

Fig. 4 .

 4 Fig. 4. Kinematics states error results.

Fig. 5 .

 5 Fig. 5. Sensors bias states error results.

  = 60 seconds period) with synthetic IMU/MAG/GNSS sensor data, the proposed CF-EKF filter is compared against the traditional centralized EKF and against the ad hoc CF/EKF filter architecture with no CF modelling compensations for a 10-run Monte Carlo simulation. Notice that the CF-EKF degenerates to the EKF when the gains k a and k m are set to zero. On the other hand, the ad hoc CF/EKF implementation is performed by setting k a and k m to zero only in the Kalman filter process matrices A and B (equations 38 and 39). This implies the same computed covariance matrices for both EKF and ad hoc implementations. Kalman filter is tuned by means of nominal sensor statistics in all three cases. Similarly, complementary gains k a and k m are kept constant throughout different implementations. Finally, the arbitrarily chosen orientation of the simulated vehicle is such that the yaw ψ, pitch θ and roll φ angles are given by

			δp(t), δv(t), ψ(t) (correction)
		ω(t)		q i c (t)
	Basic	f (t)	Complementary	p c i (t)
	Sensors		Filter		Extended
		b(t)		i v c i (t)	Kalman
					Filter
			Aiding
			Sensors
	∇(t), ε(t), ∆(t) (calibration)
	Fig. 3. CF-EKF filter overall architecture.
	white nor Gaussian in reality. Therefore, it is the only noise
	component that might require trial-and-error tuning and
	careful consideration of vehicle nominal trajectories.
		5. SIMULATION RESULTS
	In a simulated planar circular trajectory (10m radius and
	T ψ(t) = θ(t) = φ(t) =	2π T	t	(42)
	Figures 4 and 5 represents position, velocity, misalignment,
	accelerometer bias, rate-gyro drift and magnetometer bias
	root mean square (RMS) estimation errors and their
	respective predicted error covariance values in function of
	time.
					36,
					one can straightforwardly tune the EKF stage by means
					of manufacturer sensor specs and IMU sampling rate.
					However, the last noise component, namely, i a b b , is rarely

  On the other hand, EKF and ad hoc CF/EKF share the same error covariance values. However, the ad hoc CF/EKF RMS values rise above the 3σ curve yielding

				Kinematic states estimation errors
			3
			2
	δp i1 (m)	1
			0
			0	10	20	30	40	50	60
		0.8
		0.6
	δv i1 (m/s)	0.4
		0.2
			0
			0	10	20	30	40	50	60
			3
				EKF RMS		ADHOC RMS
	ψ 2 (deg)	2	CF-EKF RMS CF-EKF 3σ		EKF/ADHOC 3σ
			1
			0
			0	10	20	30	40	50	60
					t (s)

Proceedings of the 20th IFAC Symposium on Automatic Control in Aerospace,Sherbrooke, Quebec, Canada, August 21-25, 2016

More precisely, it is stable except for the unlikely case when B is parallel to g.