domain: phys spi

text:
Gilles Ferreres

Thierry Le Moing

Clément Toussaint

Simulation and comparison of control laws for the landing of an automonous parafoil under adverse conditions

1

Simulation and Comparison of Control Laws for the Landing of an autonomous Parafoil under Adverse Conditions

Gilles Ferreres1 , Thierry Le Moing2 and Clément Toussaint 3ONERA, Toulouse, France, 31055

A nonlinear simulation model of a parafoil was developed. Then, two pairs of flight control and guidance laws were designed. In the first approach, an MPC flight controller controls the yaw angle, and the guidance law computes a trajectory and the associated yaw angle reference signal to be followed. In the second simplified approach, the flight controller controls the yaw rate, and the guidance law computes the yaw rate reference value that minimizes the distance to the target. Each pair of flight control and guidance laws is tested in the non-linear closed loop simulation tool, first without and then with turbulence. The performance appears quite satisfactory, not only as regards the main criterion which is the distance to the target, but also relating to the magnitude of closed loop signals, such as yaw rates and asymmetric brake deflections, which remain within realistic levels.

I. Nomenclature

𝑉 ℎ

= parafoil horizontal airspeed 𝑉 𝑣 = parafoil vertical airspeed 𝑇 𝑎𝑝𝑝 = final approach time 𝑇 𝑡𝑢𝑟𝑛 = final turn time 𝑊 𝑥 , 𝑊 𝑦 , 𝑊 𝑧 = wind components in the inertial frame 𝑉 𝑤 = wind speed 𝜓 𝑤 = wind heading

II. Introduction

A. Context

Airdrop Aerial Delivery represents an efficient way to deliver equipment onto sites with no infrastructure, and without exposing cargo aircraft. However, landing accuracy of autonomous parafoil delivery systems remains an operational limitation when the mission requires a final dispersion smaller than 50m.

Designing a guidance and control system which ensures precision landing, especially in disturbed aerological environment, is a highly methodological challenge: First, parafoils are under-actuated systems mainly controlled with asymmetric brake deflection, symmetric deflection having a very limited effect on glide slope. Then, the low airspeed of parafoils compared to conventional aircraft makes the flight path very sensitive to changes in wind speed. Most notably during final manoeuvre, turbulence-induced disturbances significantly impact landing accuracy. Finally, airdrop flight control systems are equipped with a limited instrumentation usually excluding air data devices.

Achieving very stringent landing precision requirements for a wide range of weather conditions requires the implementation of innovative control and guidance algorithms efficiently exploiting the system dynamics. This objective has motivated the launch of a common pluri-annual research program between ONERA and DGA (the French Defense Agency).

The first step concerned the development of models and tools apt to represent adequately the behavior of payload/parafoil systems. A state-of-the-art 9-DoF (degrees of freedom) dynamics flight model was developed [START_REF] Toussaint | Flight dynamic modeling of the PBO parafoil using sparse preliminary flight test data[END_REF], which takes into account not only the relative rotations of the canopy and of the payload [START_REF] Barrows | Apparent mass of parafoil with spanwise camber[END_REF], but also added mass effects [START_REF] Lissaman | Apparent mass effects on parafoil design[END_REF] that are significant contributors to the canopy dynamics. Modellings found in the literature [START_REF] Barrows | Multibody parafoil model[END_REF] could provide an initial parameter set in order to represent the envisioned system, equipped with the PBO parachute used by the French army.

In a second step, instrumented flight tests were carried out by DGA, and the 9-DoF model was updated to match the measured trajectories and flight parameter evolutions during piloted glides. The updated model was then implemented in simulation tools aimed at assessing the factors limiting landing accuracy of autonomous parafoil landing systems, through realistic scenarios (mostly weather scenarios). Preliminary control and guidance laws were developed, and were confirmed to be critical to the achieved performance and robustness.

Additionally, with the aim of gaining experience in the field of parafoil control systems, ONERA engaged in experimental studies with a mini-paramotor. This test-bed is currently used to try and validate new landing guidance and control concepts.

The most common approaches encountered in the literature address the problem of terminal guidance in a twostep strategy:

1) A terminal path planner delivers an optimal flyable path given a known wind profile 2) The optimal trajectory is tracked by a specific controller which commands canopy brake deflections.

The trajectory planning performance is limited by the capability to predict wind along the final trajectory, and by the approximations of the model used in the path planner, which is generally a kinematic model. Moreover, the low control capacity limits the opportunity to correct for guidance errors and wind disturbances. Therefore, the optimal trajectory needs to be periodically updated during the final manoeuvre, accounting for the current state and wind estimation.

Numerous authors have proposed a wide range of optimal parafoil guidance strategies:

1) The band limited method presented in [START_REF] Carter | Band limited guidance and control of large parafoils[END_REF] encompasses the actuator limitation in the optimization problem 2) Modified Dubins paths are proposed in [START_REF] Rademacher | In-flight trajectory planning and guidance for autonomous parafoils[END_REF][START_REF] Van Der Kolf | Flight Control system for an autonomous parafoil[END_REF] for terminal guidance phase planning with a concept of altitude margin 3) Flight path parameterization using single or connected cubic Bézier curves was shown to be an attractive path-planning scheme, offering geometrically flexible trajectory shapes that maintain landing accuracy in the presence of terrain obstacles [8] 4) An optimal terminal guidance strategy based on the inverse dynamics in the virtual domain is proposed in [START_REF] Yakimenko | Using Direct Methods for terminal guidance of autonomous aerial delivery systems[END_REF][START_REF] Slegers | Optimal Control for Terminal Guidance of Autonomous Parafoils[END_REF][START_REF] Slegers | Terminal Guidance of Autonomous Parafoils in High Wind-to-Airspeed Ratios[END_REF]. This efficient algorithm which can be used in real time accounts for varying wind and disturbances by periodically updating the trajectory during the final turn. All these techniques provide a reference signal for the heading rate or heading angle. The aim of the flight controller is to generate brake deflections that track this reference signal, i.e. the true heading rate or angle should (exactly) follow this reference signal.

Among the variety of tracking controllers proposed in the literature, Model Predictive Control (MPC) seems to be an efficient way to control the yaw angle of the canopy over the terminal trajectory horizon [START_REF] Slegers | Model Predictive Control of a Parafoil and Payload System[END_REF] [START_REF] Slegers | Optimal Control for Terminal Guidance of Autonomous Parafoils[END_REF][START_REF] Slegers | Terminal Guidance of Autonomous Parafoils in High Wind-to-Airspeed Ratios[END_REF][START_REF] Alaniz | Model predictive control with application to real time hardware and a guided parafoil[END_REF]. The MPC controller uses an internal model of the turn rate dynamics to determine an optimal set of brake deflection commands, given a heading reference path.

A simple two DoF model of the roll and yaw dynamics was used to develop a yaw controller [START_REF] Slegers | Optimal Control for Terminal Guidance of Autonomous Parafoils[END_REF][START_REF] Slegers | Terminal Guidance of Autonomous Parafoils in High Wind-to-Airspeed Ratios[END_REF]. This fourth order linear model was obtained by linearizing the 6-DoF model assuming a constant air velocity.

Ward and Fowler showed that a single degree-of-freedom linear model including a turn rate first order response may be used quite effectively for parafoil inner-loop tracking [START_REF] Ward | Adaptive Glide Slope Control for Parafoil and Payload Aircraft[END_REF].

In short, the rationale for focusing on the terminal guidance phase is that the sought landing accuracy relies chiefly on its performance, in particular on its capacity to compensate efficiently for initial position errors and unknown wind perturbations until touchdown.

B. Content

This paper focuses on the terminal approach phase, which consists in a final turn. The flight controller and the guidance law are designed to ensure that the parafoil lands as close as possible to the landing target in the presence of turbulence.

Realistic simulation models and tools are a pre-requisite to convincingly address this topic. Based on previous studies [START_REF] Toussaint | Flight dynamic modeling of the PBO parafoil using sparse preliminary flight test data[END_REF] funded by DGA, in which specific instrumented piloted parafoil flight tests were performed, various models have been developed in order to gain a fine knowledge of parafoil flight dynamics, and to support the design of flight controllers. Although the reference was shown to be the 9 DOF model, simpler 6-DOF models generally turn out accurate enough to assess guidance and control laws.

The results presented in this paper concern a mini-paramotor employed in experimental studies. First, a representative 6-DOF model has been developed and integrated in a non-linear simulation tool. Second, two pairs of flight control and guidance laws have been developed. In the first approach, an MPC flight controller controls the yaw angle, and the guidance law computes a trajectory and the associated yaw angle reference signal to be followed. This guidance law makes use of wind-fixed frame [START_REF] Rademacher | In-flight trajectory planning and guidance for autonomous parafoils[END_REF], which greatly simplifies the trajectory planning formulation. In the second simplified approach, the flight controller controls the yaw rate, and the guidance law computes the yaw rate reference value that minimizes the distance to the target. Both guidance laws exploit a kinematic model, instead of a dedicated model matching exactly the dynamics of the inner closed loop corresponding to the flight controller.

Last, each pair of flight control and guidance laws is tested in the non-linear closed loop simulation tool, first without and then with turbulence. The performance appears quite satisfactory, not only as regards the main criterion which is the distance to the target, but also relating to the magnitude of closed loop signals, such as yaw rates and asymmetric brake deflections, which stay within realistic levels.

III. Experimental setup

A research testbed was developed on the basis of a mini-paramotor designed by the company Opale Paramodels. It aims at studying the dynamic behavior of a parafoil-load system, and at investigating Guidance, Navigation and Control (GNC) concepts.

It is employed in an on-going ONERA research project, whose aim is to improve the accuracy of terminal guidance of autonomous parafoil systems and to develop experimental skills.

The mini-paramotor consists of a trike, a canopy, and a brushless motor connected to a 15" propeller (only used to make the mini-motor gain altitude), and two servomotors to activate the canopy brake lines (Fig. 1).

With the canopy area of about 2.4m², the system is capable of carrying a total mass of about 6-7 kg. The forward airspeed is around 8m/s and the descent rate around 2.4 m/s. A specific avionics system has been specified and implemented for this project, which integrates off-the-shelf components, including:

1) On board computer 2) Inertial measurement unit 3) GNSS 4) Barometer 5) Anemometer 6) Servomotors 7) Potentiometers for copying positions of brake lines.

Fig. 1 Mini-paramotor and payload system.

An initial avionics had been based upon the 3DR Pixhawk, which facilitates the integration of all sensors allowing the design of a flight controller. A new type of avionics has been developed, based on a Raspberry Pi 3 and a ROS architecture.

Fig. 2 GNC hardware.

This single-board computer offers higher cpu performance and memory capacity, necessary to implement advanced path planning algorithms. On the other hand, the hardware capabilities of the Raspberry Pi's interface are not adequate for building a flight controller, so that an EMLID Navio 2 shield was used to facilitate the development of a flight controller. This shield is equipped with double IMU, Ublox GNSS, barometer, ADC interface, PWM output to control servomotors and SBUS input to decode RC signals.

Internal IMU and GNSS have been replaced by external devices for better accuracy. The avionics is currently fitted with a SBG Ellipse-A IMU and a Drotek RTK GNSS (Fig. 2).

In its standard configuration, the mini-paramotor is remotely controlled by the pilot using a radio transmitter. In this case, the radio receiver directly controls the propeller engine and the brake line servomotors.

An Sbus signal is also sent to Raspberry Pi through the Navio 2 shield, to handle flight control modes and to control data logging. An electromechanical switch, controlled by radio transmitter, allows the servomotor control from either the radio receiver, or from the flight control laws implemented in the Raspberry Pi.

The mini-paramotor is also equipped with a remote identification device, required by national regulations for all unmanned aircraft weighing more than 800g.

A ground station, based on the open source software QGroundControl running on Linux PC, provides a bidirectional radio modem link to the vehicle.

The MAVLink communication protocol is used to monitor key flight parameters, and also to send base station RTK corrections to the mini-paramotor. The ROS flight control software performs the following functions:

1) Acquisition of sensor data 2) Sbus signal decoding 3) Kinematic state estimation 4) Wind estimation 5) Terminal path planning 6) Tracking control 7) Servomotor commanding (Hitec HS-785 HB sail winches) 8) Communication with the ground station 9) Data logging to a USB memory stick.

The kinematic state estimator designed by ONERA makes use of invariant Kalman filtering. An additional estimator was developed to provide a wind velocity estimate from measurement of GNSS ground velocities and measurement of an anemometer used during the non-propelled flight phases. This configuration does not allow simultaneous estimation of the three components of the wind velocity, and remaining uncertainties in the wind estimation must be compensated by an adaptive trajectory planning.

On-line estimation of key parameters of the parafoil (air speed, descent rate, glide ratio) during the flight phases before the final turn can also improve the accuracy of final path planning.

Sensor drivers are implemented as separate ROS nodes, and a specific ROS architecture has been designed to provide a strict task scheduling. Thus, all flight control system functions are synchronized by the acquisition of IMU data at a rate of 50 Hz.

ROS messages including sensors measurements and computed data are logged to binary files on a USB memory stick. These files can be directly loaded to Matlab workspace for data analysis.

IV. Parafoil model A. Identification of the parafoil model

Most flight dynamics simulations of parafoil systems rely on 6 or 9-DoF models, depending on whether relative payload motion is considered or not. A 6-DoF model was retained for this design and for the assessment of terminal guidance and control laws, since literature and our experience show that a 6-DoF model yields good enough agreement with experimental data and is well-suited for control law design.

The structure of the 6 DoF parafoil dynamic model adopted in this work is similar to that of an aircraft. The parafoil payload system is assumed to be a rigid body. The equations of motion of the parachute-payload system are written in the body reference frame, by summing up forces and moment contributions about the system center of gravity (CG). Aerodynamic modelling consists of static and dynamic derivatives, which are constant values. Unlike in aircraft flight dynamics modelling, the air velocity considered in the aerodynamic model must be evaluated at the aerodynamic center of the canopy, and not at the center of gravity.

Aerodynamic derivatives have been estimated from flight test experiments, remotely conducted for flight data acquisition with the avionics system.

Flight tests are performed in restricted areas reserved for ONERA flight experiments. Gliding phase duration is limited due to the maximum flying altitude permitted by the UAS regulations, 120m above ground level.

In a first step, processing of steady state flights performed with different symmetric brake inputs led to estimate the parameters of lift and drag models.

In a second step, analysis of steady-state turns with asymmetric brake deflections enabled estimation of the main lateral aerodynamic derivatives. Last, a system identification technique based on an output error method was applied to update the model parameters, using flight tests performed with dynamic brake inputs.

Figure 3 presents an illustrative example of a comparison between outputs of the identified model and flight test measurements. This result highlights the limitation of the 6-DoF parafoil model: yaw rate time histories exhibit oscillations with period around 1s, due to the yawing motion of the trike supporting the avionics. The chosen model structure does not allow restitution of this yawing mode. A specific open-loop process was settled to identify a dynamic model of the servomotors installed in the miniparamotor. In order to have no unexplained and spurious delays, as it is often the case with applications such as ROS, a dedicated ROS node was developed to perform the generation of a chirp signal, the sending to PWM servo driver, the acquisition of servomotor position and the recording of all data.

It was checked that the delay introduced by the real-time execution on the Raspberry Pi 3 is around 1ms, and so negligible compared to a sampling frequency fs = 50Hz.

Tests were carried out with 4 magnitudes, PWM pulses of 50, 100, 200 and 400us respectively, to check the dependency of the results on magnitude. With a simple system like a servomotor, using identification methodologies available in toolboxes is suitable. The CAPTAIN toolbox [START_REF] Young | Recursive estimation and time-series analysis: an introduction for the student and practitioner[END_REF] was chosen to identify a continuous-time (CT) model of the servomotors, since this toolbox was utilized to identify electromechanical systems [START_REF] Young | Recursive estimation and time-series analysis: an introduction for the student and practitioner[END_REF][START_REF] Janot | Identification and control of electro-mechanical systems using state-dependent parameter estimation[END_REF][START_REF] Brunot | An improved instrumental variable method for industrial robot model identification[END_REF].

The time domain refined instrumental variable (riv) technique implemented in CAPTAIN was used to identify the structure and parameters of the following hybrid Box-Jenkins model (BJ model):

𝑦(𝑡 𝑘 ) = 𝑒 -𝑛 𝑑 .𝑇 𝑒 .𝑠 𝐵(𝑠) 𝐴(𝑠) 𝑢(𝑡 𝑘 ) + 𝐷(𝑧 -1 ) 𝐶(𝑧 -1 ) 𝜀(𝑡 𝑘 ) (1)
Where:

𝐴(𝑠) = 𝑠 𝑛 𝑎 + 𝑎 𝑛 𝑎 -1 𝑠 𝑛 𝑎 -1 + ⋯ + 𝑎 1 𝑠 + 𝑎 0 𝐵(𝑠) = 𝑏 𝑛 𝑏 𝑠 𝑛 𝑏 + 𝑏 𝑛 𝑏 -1 𝑠 𝑛 𝑏 -1 + ⋯ + 𝑏 1 𝑠 + 𝑏 0 (2)
𝐶(𝑧 -1 ) = 𝑐 𝑛 𝑐 𝑧 -𝑛 𝑐 + 𝑐 𝑛 𝑐 -1 𝑧 -𝑛 𝑐 +1 + ⋯ + 𝑐 1 𝑧 -1 + 𝑐 0 𝐷(𝑧 -1 ) = 𝑑 𝑛 𝑑 𝑧 -𝑛 𝑑 + 𝑑 𝑛 𝑑 -1 𝑧 -𝑛 𝑑 +1 + ⋯ + 𝑑 1 𝑧 -1 + 𝑑 0 𝑛 𝑎 and 𝑛 𝑏 are the degrees of the polynomials 𝐴(𝑠) and 𝐵(𝑠), 𝑛 𝑐 and 𝑛 𝑑 are the degrees of the polynomials 𝐶(𝑧 -1 ) and 𝐷(𝑧 -1 ), 𝑛 𝑑 is the number of delays and 𝑇 𝑒 is the sampling period.

Results show that a hybrid BJ model with 𝑛 𝑎 = 3, 𝑛 𝑏 = 2, 𝑛 𝑐 = 3 and 𝐷(𝑧 -1 ) = 1 gives responses in good agreement with the experiment, and the CT model is weakly magnitude dependent. Figures 4 and5 display the bode diagrams of the four CT models identified for the left and right servomotors.

To assess the quality of the simulation model, a comparison is shown in Figures 6 and7 between the time-domain model response and the measured position.

V. Path planning A. Introduction

The objective of the terminal path planning is to find the optimal turn that satisfies boundary conditions defined by the current and final states at touchdown. This phase is the last chance to manoeuvre the vehicle, and the terminal guidance scheme must compensate for initial position and other errors. The remaining flight time is unknown, and depends on initial altitude, descent rate and vertical wind velocity. Correct assumptions about the wind are crucial for the landing precision, as wind drift is the result of the integration of the wind profile along the remaining trajectory.

The determination of the trajectory is approached as a parametric optimization problem, including a kinematic model of the parafoil to reduce computational burden and to allow on-line optimization. The methodology used considers the following assumptions:

1) Initial state is defined by the cartesian coordinates (𝑥0, y0, z0) and heading 2) Landing is into the wind to reduce relative speed to ground and to improve longitudinal landing accuracy 3) The trajectory ends with a final straight line into the wind to stabilize payload motion. 4) Wind speed 𝑉 𝑤 is unknown, but assumed to be constant during this manoeuvre 5) In first approximation, horizontal airspeed 𝑉 ℎ and descent rate 𝑉 𝑣 are also assumed to be constant.

The final guidance setup depicted in Fig. 8 is similar to those considered in Yakimenko and Slegers [START_REF] Yakimenko | Using Direct Methods for terminal guidance of autonomous aerial delivery systems[END_REF]. It starts with a turn, followed by a straight section into the wind.

The final approach time 𝑇 𝑎𝑝𝑝 is chosen to stabilize the payload swinging motion and to prepare for the flare manoeuvre.

The nominal turn start position is determined on the basis of the wind velocity, assuming a constant turn rate. First, the lateral distance 𝑦0 to the target is determined by the desired constant turn rate. The turn duration is deduced as:

𝑇 𝑡𝑢𝑟𝑛 = 𝜋 𝑦0 2 /𝑉 ℎ (3)
Second, the longitudinal turn start position is determined on the basis of the drift due to wind velocity, in order to ensure that the constant turn will end at the desired point:

𝑥0 = -(𝑇 𝑡𝑢𝑟𝑛 + 𝑇 𝑎𝑝𝑝 )𝑉 𝑤 + 𝑇 𝑎𝑝𝑝 𝑉 ℎ (4)
Third, assuming a zero mean velocity of the vertical wind, the sum (𝑇 𝑎𝑝𝑝 + 𝑇 𝑡𝑢𝑟𝑛 ) of the final approach time and of the turn time determines the altitude at the start of the final turn: ℎ0 = (𝑇 𝑡𝑢𝑟𝑛 + 𝑇 𝑎𝑝𝑝 )𝑉 𝑣 [START_REF] Carter | Band limited guidance and control of large parafoils[END_REF] Fig. 8 Final turn and approach.

B. Computation of the optimal trajectory

The kinematic model equations require the knowledge of the remaining time value, which is treated as a parameter of the optimization problem.

The problem equations are formulated in a wind-fixed coordinate frame, which removes wind induced time variations and simplifies the computation of a reference course angle. During a steady turn with a constant bank angle and a constant wind speed, parafoil course rate is constant in the wind-fixed frame while it is time-varying in the inertial frame.

The origin of this wind-fixed frame coincides with the origin of the inertial frame at the initial state. The target position is shifted by the wind-induced drift during the remaining trajectory:

[ 𝑥 𝑦 ] 𝜏=1 = [ 𝑥 𝑇 + 𝑇 𝑎𝑝𝑝 𝑉 ℎ 𝑐𝑜𝑠𝜓 𝑤 -𝑇 𝑡 𝑊 𝑥 𝑦 𝑇 + 𝑇 𝑎𝑝𝑝 𝑉 ℎ 𝑠𝑖𝑛𝜓 𝑤 -𝑇 𝑡 𝑊 𝑦 ] (6)
ℎ 𝜏=1 = ℎ 𝑇 + 𝑇 𝑎𝑝𝑝 𝑉 𝑣 + 𝑇 𝑡 𝑊 𝑧 [START_REF] Van Der Kolf | Flight Control system for an autonomous parafoil[END_REF] Remember 𝑇 𝑎𝑝𝑝 is the duration of the final straight line approach, 𝑇 𝑡𝑢𝑟𝑛 is the final turn time, and:

𝑇 𝑡 = 𝑇 𝑡𝑢𝑟𝑛 + 𝑇 𝑎𝑝𝑝 (8)
The coordinates of the horizontal trajectory are represented by polynomial functions of a virtual time 𝜏:

𝑥(𝜏) = 𝑎 5 1 𝜏 5 + 𝑎 4 1 𝜏 4 + 𝑎 3 1 𝜏 3 + 𝑎 2 1 𝜏 2 + 𝑎 1 1 𝜏 1 + 𝑎 0 1 (9)
𝑦(𝜏) = 𝑎 5 2 𝜏 5 + 𝑎 4 2 𝜏 4 + 𝑎 3 2 𝜏 3 + 𝑎 2 2 𝜏 2 + 𝑎 1 2 𝜏 1 + 𝑎 0 2 (10)
The polynomial order is fixed to 5, to satisfy boundary conditions that avoid discontinuities in the bank angle. The virtual time 𝜏 varies between 0 and an unknown value, which is the second parameter of the optimization problem.

Coefficients 𝑎 𝑖 𝑗 of the above polynomials are determined by the following boundary conditions:

[ 𝑥 𝑦 ] 𝑡=0 = [ 𝑥 0 𝑥 0 ] (11)
[ 𝑥ẏ ̇]𝑡=0 = [ 𝑉 ℎ 𝑐𝑜𝑠𝜓 0 𝑉 ℎ 𝑠𝑖𝑛𝜓 0 ] (12)
[ 𝑥ẏ ̇]𝑡=𝑇 𝑡𝑢𝑟𝑛 = [ -𝑉 ℎ 𝑐𝑜𝑠𝜓 𝑤 -𝑉 ℎ 𝑠𝑖𝑛𝜓 𝑤 ] (14)
[ 𝑥ÿ ]𝑡=𝑇 𝑡𝑢𝑟𝑛 = [ 𝑉 ℎ 𝑠𝑖𝑛𝜓 𝑤 -𝑉 ℎ 𝑐𝑜𝑠𝜓 𝑤 ] ( (15)
) 16
On the other hand, first and second differentiations of the polynomial functions result in:

𝑥(𝜏) = 𝑎 5 1 𝜏 5 + 𝑎 4 1 𝜏 4 + 𝑎 3 1 𝜏 3 + 𝑎 2 1 𝜏 2 + 𝑎 1 1 𝜏 1 + 𝑎 0 1 ( 17
) 𝑑𝑥(𝜏) 𝑑𝜏 = 5𝑎 5 1 𝜏 4 + 4𝑎 4 1 𝜏 3 + 3𝑎 3 1 𝜏 2 + 2𝑎 2 1 𝜏 + 𝑎 1 1 ( 18
) 𝑑 2 𝑥(𝜏) 𝑑𝜏 2 = 20𝑎 5 1 𝜏 3 + 12𝑎 4 1 𝜏 2 + 6𝑎 3 1 𝜏 + 2𝑎 2 1 ( 19
)
Introducing these expressions at the initial and terminal points yields a system of linear equations for polynomial coefficients.

Solving this system of equations allows the computation of position and speed at each time in the virtual domain [0 𝜏 𝑡𝑢𝑟𝑛 ]. It is worth noting that this solution does not guarantee a constant speed along this path.

Then, the trajectory is sampled over a set of regularly spaced points in the virtual domain. The time interval between two consecutive points is calculated, so as to move along the trajectory with a constant speed 𝑉 ℎ :

Δ𝑡 𝑘 = 𝑑𝜏 √ 𝑑𝑥(𝑘.𝑑𝜏) 𝑑𝜏 2 + 𝑑𝑦(𝑘.𝑑𝜏) 𝑑𝜏 2 𝑉 ℎ ⁄ (20)
This time step Δ𝑡 𝑘 is used to compute the altitude change over this interval:

Δℎ 0 = 𝑉 𝑣 . Δ𝑡 𝑘 (21)
The reference course angle is computed from the horizontal speed components:

ψ 𝑘 = atan ( 𝑑𝑦(𝑘.𝑑𝜏) 𝑑𝜏 𝑑𝑥(𝑘.𝑑𝜏) 𝑑𝜏 ⁄ ) ( 22
)
and the course rate is simply evaluated by numeric differentiation:

𝑑𝜓 𝑘 𝑑𝑡 = (ψ 𝑘 -ψ 𝑘-1 )/Δ𝑡 𝑘 ( 23
)
These values are used to evaluate a cumulative performance index that must be zero for the optimal trajectory:

𝐶 = 𝜇 1 (∑ 𝑉 𝑣 𝑘 . Δ𝑡 𝑘 -ℎ 0 ) 2 + 𝜇 2 (24)
This quadratic criterion can be completed by an additional constraint allowing the limitation of the course rate: Note that the parachute dynamics cannot be taken into account with this calculation scheme. However, it is possible to introduce quasi-static phenomena, such as the rate speed change observed in steady turn, which can be modelled by the equation: The optimization problem is solved using a Gauss-Newton algorithm, which generally converges after a few iterations as long as a rigorous solution exists. During the periodic update of the optimal trajectory, the initial value of the parameters is set from the solution of the previous optimization, which speeds up convergence.

𝑉 𝑣 (𝑘) = 𝑉 𝑣 + 0.25 ( 𝑑𝜓 𝑘 𝑑𝑡 ) 2 (21)

C. Some results

Figure 10 shows an example of the optimal final turn obtained with a wind of 5m/s. This result exhibits the benefit of the path optimization in the wind fixed frame. The turn rate of the trajectory relative to the wind (in blue) indicates two maxima close to 25°/s, whereas the ground turn rate (in red) is much higher during the last quarter turn.

Fig. 10 Example of optimal flight path in air frame and inertial frame.

The path planning algorithm can generate reference trajectories for a range of initial conditions, as shown in Fig. 11. Ten trajectories are computed with a random selection of initial position (𝜎 𝑥 =𝜎 𝑦 =7m, 𝜎 𝑧 =2m) and heading (𝜎  =10°). It can be observed that these trajectories still converge to the target.

When the altitude is a little too low with respect to the horizontal distance, the optimal trajectory requires turn rates reaching 40°/s. MPC control could have difficulties to produce an accurate tracking for these high turn rates, so that it may be necessary to optimize the trajectory with a turn rate boundary. Figure 12 shows new trajectories generated with the same initial points, but with the introduction of a turn rate limit of 30°/s.

In some cases, no optimal trajectory reaches the target. As depicted by the first plot, the landing position error still remains on the longitudinal axis.

VI. Model Predictive Control

A. Principle Model predictive control (MPC) is a well-established approach for advanced system control in many industrial applications. MPC is efficient, simple in tuning and robust for tracking a known reference trajectory considering constraints and disturbances. This method is often used in parafoil guidance, navigation and control systems, to track a reference yaw or course angle provided by the path planning. It was successfully demonstrated by Slegers and Costello [START_REF] Slegers | Model Predictive Control of a Parafoil and Payload System[END_REF].

The MPC methodology attempts to solve an on-line open loop finite horizon optimal control problem subject to input, state, and/or output constraints. As shown in Fig. 13, at time t, the system model is used to predict the future behavior of the controlled plant over the prediction horizon Np.

Fig. 13 Example of optimal flight path in air frame and inertial frame.

The method delivers an optimal control input signal u(t) over a control horizon Nc (usually Nc ≤ Np), which minimizes the difference between the predicted output and the reference trajectory. Only the first element of the optimal signal is applied to the plant, and all other elements are discarded. Then, at the next sampling instance, the whole procedure is repeated.

MPC allows operation within constraints, and it can be extended to multi-variable and nonlinear systems. However, the method has some limitations. As the system complexity increases or with non-linear models, the on-line calculation burden is substantially increased. It could also be difficult to predict the controller performance and robustness with the real system in closed loop.

With a linear model and without constraints, the control law can be computed analytically as the solution of a cost optimization problem. Consider the following discrete-time linear system in state space form:

𝑥 𝑘+1 = 𝐹𝑥 𝑘 + 𝐺𝑢 𝑘 𝑦 𝑘 = 𝐶𝑥 𝑘 ( 26
)
The finite horizon optimal control problem to be solved at each sampling instance can be expressed as the following unconstrained optimization problem. The optimization criterion J is defined as:

𝐽 𝑘 = (𝑌 𝑘+1 𝑟 -𝑌 𝑘+1 ) 𝑡 𝑄(𝑌 𝑘+1 𝑟 -𝑌 𝑘+1 ) + 𝑈 𝑘 𝑡 𝑅𝑈 𝑘 ( 27
)
The input signal to be optimized is:

𝑈 𝑘 = [𝑢 𝑘 , 𝑢 𝑘+1 , … , 𝑢 𝑘+𝐻-1 ] 𝑡 ( 28
)
where H is the prediction horizon. The obtained output signal is:

𝑌 𝑘+1 = [𝑦 𝑘+1 , 𝑦 𝑘+2 , … , 𝑦 𝑘+𝐻 ] 𝑡 ( 29
)
while the reference trajectory is:

𝑌 𝑘+1 𝑟 = [𝑦 𝑘+1 𝑟 , 𝑦 𝑘+2 𝑟 , … , 𝑦 𝑘+𝐻 𝑟 ] 𝑡 (29)
Matrices Q and R are used to balance the respective minimization of the output error and control input costs. The discrete-time linear model can be used to estimate the future output over the prediction horizon:

𝑌 𝑘+1 = 𝐾 𝐶𝐹 𝑥 𝑘 + 𝐾 𝐶𝐹𝐺 𝑈 𝑘 ( 30
)
where:

𝐾 𝐶𝐹 = [ 𝐶𝐹 𝐶𝐹 2 ⋮ 𝐶𝐹 𝐻 ] ( 31
) 𝐾 𝐶𝐹𝐺 = [ 𝐶𝐺 𝐶𝐹𝐺 𝐶𝐹 2 𝐺 ⋮ 𝐶𝐹 𝐻-1 𝐺 0 𝐶𝐺 𝐶𝐹𝐺 ⋮ 𝐶𝐹 𝐻-2 𝐺 … … … ⋱ … 0 0 0 ⋮ 𝐶𝐹𝐺 0 0 0 ⋮ 𝐶𝐺 ] (32)
The optimal control signal 𝑈 𝑘 over the prediction horizon, which minimizes the cost 𝐽 𝑘 , is:

𝑈 𝑘 = (𝐾 𝐶𝐹𝐺 𝑡 𝑄𝐾 𝐶𝐹𝐺 + 𝑅) -1 𝐾 𝐶𝐹𝐺 𝑡 𝑄(𝑌 𝑘+1 𝑟 -𝐾 𝐶𝐹 𝑥 𝑘 ) (22)

B. Application to the parafoil

The prediction model used in the above MPC algorithm was identified from the yaw response of the 6-DoF parafoil model. Fig. 14 shows that a good approximation can be obtained with a reduced third order model in the following form:

𝑥̇= [ 0 1 0 0 0 1 0 -𝜔² -2𝜔𝜁 ] 𝑥 + [ 0 𝑏 𝑎 𝜔² ] 𝑢 (33) 𝜓 = [1 0 0]𝑥 ( 34
)
where 𝑎 = -1.15, 𝑏 = 0.2, 𝜔 = 1.33 𝑟𝑎𝑑/𝑠 and 𝜁 = 0.32. This continuous-time model is discretized inside the MPC algorithm.

The output of the prediction model is quite close to the yaw output (3 rd subfigure of Fig. 14, from above) and to the course angle (4 th subfigure) of the 6-Dof model.

During turn, these flight parameters are not identical even without wind. This can be illustrated by the result of Fig. 15, obtained by applying the MPC algorithm to the 6-DoF model, to track a reference output trajectory generated by the path planning algorithm.

The difference between yaw and course angles reaches more than 3°. It is mainly due to the effect of angle of attack when the parafoil is banked, while the angle of side slip remains low. Last note that the MPC algorithm can be quite easily extended to deal with limitations on brake line commands or parafoil turn rate. A quadratic programming approach can be used with these linear inequality constraints. Penalty functions can be incorporated into the optimization criterion to simplify the implementation. Figure 16 illustrates the tracking of the course angle after the introduction of a turn rate limit set at 20°/s.

VII. Simulation results

Performance of the proposed path-planning and control algorithms is demonstrated through detailed example simulations and Monte Carlo studies in order to evaluate the effects of turbulence.

A. Context

Common turbulence models used in simulation scenario model are Dryden and Von Karman models. The continuous Dryden model was chosen because of its easier implementation. According to the military references [START_REF]US Military Specification: Flying Qualities of Piloted Airplanes[END_REF], the turbulence scale lengths and the turbulence intensities are functions of the altitude. The standard deviation of the wind was directly imposed in this evaluation. To generate a signal with correct characteristics, a band-limited white noise signal is passed through low-pass filters, derived from the spectral square roots of the Dryden power spectrum functions.

Flight path generation, MPC control and non-linear parafoil flight dynamics model have been implemented in a comprehensive simulation framework for algorithm validation and performance evaluation.

Simulation of energy management and homing phases are not implemented, but it is assumed that they are performed imperfectly (due to wind shifts or other factors), and thus initialization of terminal guidance may not be performed in an ideal geometry. In this case, the parafoil begins terminal guidance heading downwind, where the wind blows to the north at 2 m/s.

Path planning requires wind knowledge along the remaining trajectory. The deterministic wind component was assumed constant and the estimation algorithm was not modeled: the wind velocity used by trajectory optimization was delivered by low pass filtering of the instantaneous wind speed.

B. A result without turbulence

Figure 17 shows an example of a reference path tracking with a constant wind of 2m/s and without turbulence. The parafoil initially tracks the optimal trajectory. However, a deviation of a few meters appears after a couple of seconds and the impact point is 3m away from the target. This final error can be explained by the air velocity variations due to the low damping of the phugoid mode, as depicted in Fig. 18. The average value of the oscillations is also different from the equilibrium value used to compute the reference trajectory.

The implementation of an output feedback stabilization law which acts on symmetric deflection of the canopy trailing-edge strongly reduces airspeed oscillations and thus the impact error, as shown in Fig. 19 and20.

C. Results in the presence of turbulence

To illustrate the effect of wind variation, Fig. 21 shows the results of one simulation performed with turbulence. The Dryden turbulence model described above was used with a horizontal wind standard deviation of 0.5 m/s. Impact error is greater than 10 meters.

This landing error can be reduced by updating the reference trajectory en route. Fig. 22 shows the new trajectory obtained with the same turbulence and a replanning each second, enabling the parafoil to land 3m from the target. A series of Monte Carlo simulations were performed to assess the overall effect of turbulence. Figure 23 presents the landing dispersion obtained without trajectory replanning. Longitudinal impact errors vary between -15m and +15m, and lateral impact errors are approximately two times lower.

Dispersion can be characterized by the elliptical error probability (EEP), which defines the radius of the ellipse including a given percentage of the impacts. EEP encompasses a more general definition of circular error probability (CEP) when variance is unequal in the dimensional axes.

Dispersion results are shown in Fig. 23. The resulting circular error probability (CEP) shown by the circle is 16.8m and is defined by the radius which includes 50 per cent of the impacts.

Activation of reference trajectory updating each second reduces the ellipse radii by a factor of about 3, as shown in Fig. 24.

Figure 25 shows results of a new Monte Carlo simulation, carried out without the longitudinal control law. It generates an impact offset of 4m in the opposite direction of the wind, and landing dispersion increases by approximately 25%.

VIII. A simplified guidance and heading rate control approach

In this section the flight controller controls the yaw rate, instead of the yaw angle above, and the guidance law computes the yaw rate reference value that minimizes the distance to the target on the kinematic model, instead of a yaw angle reference signal above. As this yaw rate reference value is constant, although it may be re-adjusted each second, the flight controller need not be predictive, its aim is only to impose a constant value for the yaw rate output.

Figure 26 illustrates the principle of this method. The position of the starting point is defined to have an initial trajectory (in dotted line) passing through this point and the target point. Trajectory updating provides a new circular arc (in solid line) starting at the current position, its radius is chosen to minimize the distance to the target.

A. Longitudinal airspeed control

As a preliminary, a longitudinal flight control law is synthesized, which uses the symmetric brake deflection input to control the airspeed 𝑉. When using the kinematic model for the design of the guidance law below, the horizontal and vertical speeds are assumed to be constant, which is not true in practice even without turbulence and when maintaining the symmetric brake deflection input to a constant value. Indeed, a high value of the asymmetric brake deflection input (used by the lateral flight control law) induces a strong coupling on the longitudinal variables, and thus a significant variation of the speeds. Since it is not possible to regulate both the horizontal and vertical speeds, a realistic choice is to regulate the airspeed.

Since there are 4 states in the longitudinal model, obtained by linearizing the nonlinear 6 DOF model of the parafoil at a trim point, 4 longitudinal outputs 𝑎𝑧 (vertical acceleration), 𝑞 (pitch rate), 𝜃 (pitch angle), 𝑉 and an integrator on 𝑉 are used to place the 5 main closed loop poles, namely the integrator pole and the short period and phugoid modes. The unplaced actuator closed loop poles are checked a posteriori. Classically, the frequencies of the closed loop short period and phugoid modes are chosen to be the same as the open loop ones, and the damping ratio is fixed to 0.7. One does not accelerate the longitudinal closed loop, to reduce the magnitude of the symmetric brake deflection input signal, since there is a limited authority on this signal.

B. Lateral heading rate control

The architecture of the flight controller is:

𝑢 = -𝐾[∫(𝜓 ̇-𝜓 ̇𝑟𝑒𝑓 ) 𝜓 ̇𝑝] + 𝐻𝜓 ̇𝑟𝑒𝑓 ( 35
)
Where u is the plant input (asymmetric brake deflection), 𝜓 ̇ is the yaw rate, 𝜓 ̇𝑟𝑒𝑓 is the yaw rate reference value and p is the roll rate. K is the feedback gain, and H a static feedforward term. K need not be computed here, since it can be extracted from the MPC flight controller above, which contains a feedback part. Otherwise, different classical methods could be used to design this state or output feedback controller, e.g. a modal technique which places the main closed loop poles on the basis of the open loop plant model (see the previous subsection). This one can be the 3rd order simplified model used by the MPC flight controller, or a higher order model obtained by linearizing the nonlinear 6 DOF model of the parafoil at a trim point.

The feedforward term H is a tuning parameter that will be chosen in the following to minimize the distance to the target.

C. Computation of the heading rate reference value

The equations of the kinematic model are recalled first:

𝑥̇= 𝑉 ℎ cos(𝜓) + 𝑉 𝑊 cos (𝜓 𝑊 )

𝑦̇= 𝑉 ℎ sin(𝜓) + 𝑉 𝑊 sin (𝜓 𝑊 )

𝑧̇= 𝑉 𝑣 (38)

The horizontal and vertical speeds 𝑉 ℎ and 𝑉 𝑣 are assumed to be known and constant, as well the constant wind defined by 𝑉 𝑊 and 𝜓 𝑊 . Let 𝑥 0 ,𝑦 0 ,𝑧 0 ,𝜓 0 be the initial values of x,y,z,𝜓. If 𝜓 = 𝜔𝑡 + 𝜓 0 , the solution to the differential equations above is: The issue is to minimize the distance to the target with respect to 𝜔 between -𝜔 𝑚𝑎𝑥 and 𝜔 𝑚𝑎𝑥 , where 𝜔 𝑚𝑎𝑥 is a tuning parameter. This approach is suboptimal in the sense that it may not be possible to minimize this distance to zero.

𝑥(𝑡) = 𝑥 0 + 𝑉 ℎ 𝜔 (sin(𝜔𝑡 + 𝜓 0 ) -sin (𝜓 0 )) + 𝑉 𝑊 cos(𝜓 𝑊 ) 𝑡 (
Remark: there is no constraint on the final value of the heading angle in this simplified approach, whereas this final value must be 𝜓 𝑊 + 𝜋 in the first approach.

D. Comparison of simulation results with the first and simplified approaches

A preliminary step is to tune the value of the static feedforward term H. To this aim, the distance to the target is calculated, as a function of H, on the nonlinear simulator without turbulence. The following table gives the results, where the first line is the value of H, the second line is the distance to the target, and the third line is the normalized magnitude of the input signal u (asymmetric brake deflection), constrained between -1 and 1: With H=-3, in the presence of turbulence, the average value of the distance to the target is 4.77 m and the standard deviation is 3.20 m. A better result is obtained with the first approach: the distance to the target is 0.17 m without turbulence. In the presence of turbulence, the average value is 3.24 m and the standard deviation is 2.68 m. Nevertheless, the results of the simplified approach are not so far.

As for the magnitude of the heading rate in the presence of turbulence, its average value (standard deviation) is 31.75 deg/s (3.84 deg/s) with the first approach and 33.60 deg/s (3.27 deg/s) with the second one. The results are close, and reasonable.

Concerning the magnitude of the input signal u in the presence of turbulence, its average value (standard deviation) is 0.61 (0.04) with the first approach. But the values are significantly higher with the second approach, the average value (standard deviation) being 0.89 (0.04), which means a higher solicitation.

Last, the final value of the bank angle is studied. In the presence of turbulence, its average value (standard deviation) is -2.35 deg (6.11 deg) with the first approach and 1.52 deg (5.57 deg) with the second one. Thus, the result is slightly better with the simplified approach.

Figures. 27 and 28 (resp. 29 and 30) display the Monte Carlo landing dispersion (resp. the trajectories) with the two approaches in the presence of turbulence.

IX. Conclusion

After settling the context, the paper has presented a thorough investigation into flight control algorithms for parafoil automatic precision landing in turbulent weather. The flight models were based on previous studies and experiments involving the DGA, and then tuned to ONERA's experimental vehicle: a paramotor equipped with suitable avionics, employed as a test-bed for control laws.

Two approaches for terminal guidance and control were presented in this paper, and tested on a realistic nonlinear simulator. Satisfactory results were obtained, without turbulence and in the presence of turbulence.

The reference approach presented first relies on regularly updated optimized trajectories followed by an MPC flight controller. The trajectories have a polynomial parametrization simple enough for real-time implementation, with degrees of freedom allowing to enforce desirable landing conditions (straight flight into the wind). The trajectory is provided to the flight controller as a series of yaw angles. This approach demonstrates very small landing position errors with constant wind and varying initial positions. On a realistic scenario with turbulence, the mean errors remain well below the 50m mark typical of operational air-delivery systems. Longitudinal control on top of the MPC controller was also shown to provide enhanced accuracy thanks to smoother vertical trajectories (phugoid damping).

The second approach is a simplified algorithm based on the yaw rate control. Though it is theoretically less performing than the reference, both for following a preset trajectory and for fixing landing conditions, it demonstrated almost similar performance in terms of final accuracy. This good result can be partly explained by the fact that the simplified approach is more respectful of the parafoil dynamics since its aim is ideally to make a single turn. Being also easier to implement, this suboptimal approach would be worth further styduing.

New approaches are also being developed to provide a reference trajectory that respects the parafoil dynamics at the beginning and the end of the final manoeuvre while minimizing yaw rate variations.

All these simulation studies must be followed by the implementation of these guidance and control laws in the paramotor avionics which was designed to perform on-line optimization. Hardware-in-the-loop simulations will be conducted to validate these real-time algorithms before flight tests.

Fig. 3

 3 Fig. 3 Responses of identified lateral model.

Fig. 4 Fig. 6

 46 Fig. 4 Bode diagram of the left models. Fig. 5 Bode diagram of the right models.

Fig. 9

 9 Fig. 9 descent rate as function of turn rate.

Fig. 11 Fig. 12

 1112 Fig. 11 Optimal trajectories with initial condition errors. Fig. 12 Optimal trajectories with initial condition errors and turn rate limitation.

Fig. 14

 14 Fig. 14 Responses of prediction model. Fig. 15 Comparison between heading and course angles.

Fig. 16

 16 Fig. 16 Response of 6-DoF model with MPC and turn rate limitation.

Fig. 17 Fig. 18

 1718 Fig. 17 Simulated approach in ideal conditions (Wind 2 m/s).

Fig. 20

 20 Fig. 20 Air velocity.

Fig. 21

 21 Fig. 21 Example of simulated approach with turbulence (turb=0.5 m/s).

Fig. 23

 23 Fig. 23 Monte Carlo landing dispersion without trajectory updating.

Fig. 26

 26 Fig. 26 Geometry of the reference trajectory.

Fig. 27

 27 Fig. 27 Monte Carlo landing dispersionwith the first approach.

Fig. 28

 28 Fig. 28 Monte Carlo landing dispersionwith the second simplified approach.

Fig. 29

 29 Fig. 29 trajectories with the first approach. Fig. 30 Trajectories with the second simplified approach.

  𝑇 𝑎𝑝𝑝 𝑉 ℎ 𝑐𝑜𝑠𝜓 𝑤 -𝑇 𝑡 𝑊 𝑥 𝑦 𝑇 + 𝑇 𝑎𝑝𝑝 𝑉 ℎ 𝑠𝑖𝑛𝜓 𝑤 -𝑇 𝑡 𝑊 𝑦 ]

		[	𝑥ÿ ]𝑡=0	= [ -𝑉 ℎ 𝑠𝑖𝑛𝜓 0 𝑉 ℎ 𝑐𝑜𝑠𝜓 0	] 𝜓 0 ̇(13)
	[ 𝑦 ] 𝑥	𝑡=𝑇 𝑡𝑢𝑟𝑛	= [	𝑥 𝑇 +

  𝑇 1 = 𝑇 -𝑇 𝑎𝑝𝑝 , one chooses 𝜓 = 𝜔𝑡 + 𝜓 0 , where 𝑇 𝑎𝑝𝑝 and 𝜔 are fixed. For 𝑇 1 < 𝑡 ≤ 𝑇, the heading rate reference value is fixed to zero, to minimize the final (absolute) value of the bank angle at t = T. Thus, the distance to the target is √𝑥 2 (𝑇) + 𝑦 2 (𝑇), where: 𝜔𝑇 1 + 𝜓 0 ) -sin (𝜓 0 )) + 𝑉 𝑊 cos(𝜓 𝑊 ) 𝑇 1 + 𝑇 𝑎𝑝𝑝 (𝑉 ℎ cos(𝜔𝑇 1 + 𝜓 0 ) + 𝑉 𝑊 cos (𝜓 𝑊 )) 𝜔𝑇 1 + 𝜓 0 ) -cos (𝜓 0 )) + 𝑉 𝑊 sin(𝜓 𝑊 ) 𝑇 1 + 𝑇 𝑎𝑝𝑝 (𝑉 ℎ sin(𝜔𝑇 1 + 𝜓 0 ) + 𝑉 𝑊 sin (𝜓 𝑊 ))

				39)
	𝑦(𝑡) = 𝑦 0 -	𝑉 ℎ 𝜔	(𝑐𝑜𝑠(𝜔𝑡 + 𝜓 0 ) -𝑐𝑜𝑠 (𝜓 0 )) + 𝑉 𝑊 𝑠𝑖𝑛(𝜓 𝑊 ) 𝑡	(40)
	Let 𝑇 = . For 𝑡 ≤ 𝑥(𝑇) = 𝑥 0 + 𝑧 0 𝑉 𝑣 𝑉 ℎ 𝜔 (sin(𝑦(𝑇) = 𝑦 0 -𝑉 ℎ 𝜔 (cos(

Table . 13 Distance to the target as a function of feedforward gain H.

 .

	H=1	H=0	H=-1	H=-2	H=-3	H=-4
	22.80 m	16.02 m	10.23 m	5.27 m	0.94 m	3.19 m
	0.34	0.27	0.36	0.61	0.87	1.00

Research engineer, Information Processing and Systems Department.

Research engineer, Information Processing and Systems Department.

Research engineer, Information Processing and Systems Department.


Ali Amroune

Gilbert Fantozzi

Synthesis of Al203-SiC from kyanite precursor

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

I. INTRODUCTION

Due to their highly ordered structure, chemical stability, and high strength approaching the theoretical value, ceramic whiskers have been attracting interest for use as reinforcements in ceramic matrices for thermostructural applications. Among these, silicon carbide whiskers (SiC W ) have been extensively studied as reinforcing materials in ceramic matrices. [1][2][3][4] Three main methods have been utilized to synthesize SiC W : (i) carbothermal reduction of silica in the presence of metal catalysts including rice-hulls as starting material, [5][6][7][8][START_REF] Shalek | Proc. Int. Conf. On Whisker-and-Fiber Toughened Ceramics[END_REF][START_REF] Wang | Ceramic Powder Science[END_REF][START_REF] Shimada | [END_REF] (ii) thermal decomposition of organic silicon compounds, 12,13 and (iii) reaction between halides in presence of hydrogen. 14,15 Because of its simplicity, the first method is the most utilized. However, the toxicity (mainly by inhalation) related to whisker handling and the high cost of their production have limited their application. To reduce these disadvantages, an in situ SiC W synthesis approach has been attempted. Partial carbothermal reduction of Si 3 N 4 has been utilized 16,17 to obtain Si 3 N 4 /SiC W powders. Synthetic or natural aluminosilicates nAl 2 O 3 и mSiO 2 may give Al 2 O 3 and SiC W as final products by carbothermal reduction at high temperatures in an inert atmosphere. 18,19 It has been reported 18,19 that the fracture toughness K IC values of Al 2 O 3 -SiC W composites prepared by hot pressing from in situ synthesized powders range between 5 and 7.6 MPa m 1/2 . Similar composites prepared by conventional physical mixing of the components and sintering generally present a fracture toughness ranging between 6 and 9.5 MPa m 1/2 . These values indicate the possibility of obtaining in situ composites having mechanical properties comparable to those of similar conventional composites prepared from separately synthesized SiC W and Al 2 O 3 . In this study, carbothermal reduction of kyanite, a large availability aluminosilicate mineral with high alumina content, was used as a way to in situ obtain ␤-SiC whisker and ␣-Al 2 O 3 particles. This route would simplify the conventional costly way to prepare powder for elaborating Al 2 O 3 -SiC W composites. Some reaction conditions were studied to synthesize a powder with a favorable morphology for elaborating such composites.

II. EXPERIMENTAL

Table I lists chemical compositions of the utilized kyanite (Virginia, kyanite supplied by SOGEMET, Paris, France) as aluminosilicate precursor and indicates the content of impurities in the mineral, primarily iron and titanium oxides. X-ray analysis showed the good crystallinity of this mineral. Scanning electron microscopy (SEM) observation of the powder reveals that the particle size ranges between 5 and 53 m. As carbon sources, two types of carbon black have been utilized (supplied by DEGUSSA, Paris, France) with specific surface areas of 330 and 996 m 2 /g. XRD analysis of these two types of carbon reveals their amorphous states. Kyanite powder was mixed (in ethylic alcohol) with carbon black in a molar ratio of C/SiO 2 ≅ 5.5 by ball milling. The slurries were then dried and the powder mixtures were sieved to less 500 m.

A vertical tubular furnace (Pyrox, Paris, France) at 1260-1550 °C was used for carbothermal reduction of kyanite. A bed of the mixed starting powder (10 mm in height and 50 mm in diameter) was prepared in a graphite crucible and heat treated for2ha teach reaction temperature under argon. The heating and cooling rates were both 300 °C / h. After each heat treatment, unreacted carbon was removed by oxidation in air at 625 °C for 5 h before characterization.

X-ray diffraction (XRD; Cu K ␣ radiation with a graphite monochromator) was used for phase analysis. XRD quantitative analysis was been performed on the synthesized powders, too. Mixtures of Al 2 O 3 and SiC-whiskers in different proportions have been prepared and a calibration curve established. The morphology was investigated by SEM (JEOL 840, Tokyo, Japan) with an energy dispersive x-ray spectroscopy (EDS) system (EDAX TRACOR, Paris, France).

III. RESULTS AND DISCUSSION

A. Phase evolution XRD phase analysis (Fig. 1) shows that the carbothermal reaction goes through many steps leading to ␤-SiC and ␣-Al 2 O 3 as final products. Mullite and a small amount of an amorphous phase were detected as intermediate solid phases. The formation of ␤-SiC and ␣-Al 2 O 3 began at about 1320 and 1370 °C, respectively. Figure 1 shows that the rate of reduction of kyanite started very slowly at 1320 °C and increased with increasing temperature.

On the basis of XRD analysis the following reaction scheme is proposed. The thermal decomposition of kyanite to mullite and amorphous silica occurs below 1370 °C.

3͑Al 2 O 3 и SiO 2 ͒ kyanite → Ͻ1370 °C 3Al 2 O 3 и 2SiO 2 + SiO 2 , (1)
The silica produced by this decomposition is reduced to SiC above 1320 °C.

SiO 2 + 3C → Ͼ1320 °C SiC + 2CO , (2)
For the last reaction, several mechanisms have been proposed, all of which involve two important steps 5,7,8,20,21 The first step is related to the formation of silicon monoxide by the reaction:

SiO 2 + C → SiO + CO . ( 2a
)
This reaction is followed by

SiO + 2C → SiC + CO . (2b)
Combination of reactions (2a) and (2b) gives reaction (2). Above 1370 °C (after consumption of silica), the mullite is reduced. This reduction forms ␣-alumina (corundum).

3Al 2 O 3 и 2SiO 2 + 6C mullite → Ͼ1370 °C 3Al 2 O 3 + 2SiC + 4CO . ( 3
)
The combination of reactions ( 2) and (3) gives

3Al 2 O 3 и 2SiO 2 + SiO 2 + 9C → 3Al 2 O 3 + 3SiC + 6CO .
The reaction in the system is complete at 1550 °C.

B. Morphology of the ␣-Al 2 O 3 /␤-SiC powder

SEM observations of the powder synthesized at 1550 °C revealed that whiskers coexist with spherical particles. The synthesized powder has a woolly aspect. . Whiskers from the upper layer have a higher aspect ratio. This is likely due to the lower supersaturation degree of silicon and carbon monoxides at the upper layer. A higher supersaturation degree is not favorable to SiC whisker growth. [START_REF] Shalek | Proc. Int. Conf. On Whisker-and-Fiber Toughened Ceramics[END_REF] The mean dimensions of SiC whiskers and Al 2 O 3 particles in different parts of the reaction zone are presented in Table II.

Globules were observed at the tip of some whiskers (Fig. 3) and EDS analysis reveals that these globules are iron-rich regions with the presence of silicon. As in the case of separate synthesis of SiC whiskers, these globules are solidified droplets which form at high temperatures during the vapor-liquid-solid (VLS) 6,7 growth mechanism of SiC whiskers. The impurities, such as transition metals, constitute energetically favourable sites for SiC nucleation. The presence of whiskers having a conical tip instead of balls indicates that there was another mechanism of crystal growth. Although the nature of this mechanism is uncertain, vapor-solid (VS) growth mechanism of whiskers has been proposed, [START_REF] Campbell | Whisker Technology[END_REF] and carbon black may constitute nucleation sites for SiC. [START_REF] Kennedy | Proc. of the British Ceram. Society[END_REF] In the VS mechanism, SiC crystallisation occurs directly from the vapor phases SiO and CO formed at high temperatures.

The tip growth is postulated to occur on an axial screw dislocation, but such dislocations are rarely seen. [START_REF] Hurley | Advanced Composites Proceedings[END_REF] XRD quantitative analysis of the phases present in the powder after completion of the synthesis reaction gave approximately 28 wt% SiC in Al 2 O 3 /SiC powder. Optimal strength values for densified Al 2 O 3 /SiC W composites prepared by conventional mixing of alumina and SiC whiskers are obtained with 25 wt% SiC. [START_REF] Lio | [END_REF] This value may be achieved by adding alumina after synthesis. Adding alumina before synthesis to the raw materials does not contribute to a good morphology of the whiskers. 19 Some too-long, bent, or twisted whiskers and agglomerates are also encountered in the synthesized powders. These characteristics are also observed in the case of separate synthesized SiC whiskers. 26,27 This needs further operation of milling-deagglomeration (for example by ball milling) before sintering the powder to form composite materials with good dispersion of the phases. 28

C. Effect of carbon specific surface area

Figure 4 shows XRD patterns of heat-treated mixtures of kyanite and carbon black at the same temperature 1500 °C with two different values of Brunauer-Emmett-Teller (BET) specific surface area of carbon: 330 and 996 m 2 /g. Higher specific surface area of carbon lead to completion of the carbothermal reduction reaction, whereas a significant amount of mullite is detected in the sample with lower specific surface area at the same temperature of 1500 °C (the reaction in this latter case is completed between 1500 and 1550 °C).

Figure 5 shows the evolution of the most intense peak (111) of ␤-SiC formed from carbothermal reduction of kyanite in the temperature range 1370-1550 °C. A higher specific surface area of carbon black leads to an increase in the conversion rate to SiC. Quantitative analysis shows that the amount of SiC phase in Al 2 O 3 / SiC powder (about 28 wt%) does not change with increasing the temperature from 1500 to 1550 °C in the case of the higher surface area. A similar effect of specific surface area of carbon has been observed in the case of the carbonitriding of silica. 29,30 Higher specific surface area of carbon does not influence favorably the morphology of the synthesized powder. The SEM micrograph (Fig. 6) shows that the synthesized powder is mainly composed of agglomerates, and the fibrous phase SiC is not morphologically regular. Whiskers are bent, twisted, and very thin. This morphology is not favorable for Al 2 O 3 -SiC W composites, since these characteristics lead to large flaws, porosity, and inhomogeneous microstructure in the final composite. A carbon black with a relatively high surface area would lead to a large number of sites for the nucleation and growth of SiC during synthesis. These conditions favor irregular morphology of the whiskers.

IV. CONCLUSION

Carbothermal reduction of kyanite leads to in situ Al 2 O 3 particles and SiC-whisker formation with 28 wt% SiC. The completion of the reaction is obtained at 1550 °C when using a carbon with a specific surface area 330 m 2 /g and a heating time of 2 h. With a higher value of 996 m 2 /g, the completion of the reaction occurred at 1500 °C. During carbothermal reduction of kyanite, this mineral decomposes to free silica and mullite as intermediate phases. The reduction of mullite occurs at higher temperature than that of free silica. Mullite reduction does not occur before full consumption of the free silica. The aspect ratio of the synthesized SiC in the upper layer of the samples is relatively higher than that in the lower one. Both VLS and VS crystal growth mechanisms for SiC W are suggested. The morphology of the synthesized powder is not with a too-high specific surface area of carbon. The principal advantages of this in situ route are reducing the health hazard related to the handling of whiskers, the use of raw materials with large availability such as kyanite, and the simplification of the process of Al 2 O 3 /SiC W composites. The principal disadvantage of this method is the difficulty of morphological control of the phases because they are not synthesized separately.

FIG. 1 .

 1 FIG. 1. X-ray diffraction patterns on heat treated powders from kyanite-carbon system in the temperature range 1260 -1550 °C under argon. (A) ␣-Al 2 O 3 ; (B) ␤-SiC; (K) kyanite; (M) mullite.

FIG. 5 .

 5 FIG.5. Evolution of the most intense peak (111) of ␤-SiC formed from carbothermal reduction of kyanite in the temperature range 1370-1550 °C. With the higher value of specific surface area of carbon black the conversion to SiC is accelerated.

TABLE I .

 I Chemical composition of kyanite

	Oxydes	wt %
	SiO 2	37.64-43.70
	Al 2 O 3	54.00-60.06
	Fe 2 O 3	0.40-1.16
	TiO 2	0.52-1.65
	CaO	0.03
	K 2 O	0.42
	Na 2 O	0.42
	Weight loss	0.21

TABLE II .

 II Dimensions (m) of the synthesised SiC-whiskers and Al 2 O 3 -particles at 1550°C under argon atmosphere.

				FIG. 3. SEM micrograph showing a globule at the tip of a whisker.
				The globule is a solidified droplet which forms at high temperatures
				during the VLS growth mechanism of the SiC-whisker.
	FIG. 2. SEM micrographs showing difference in morphology between
	reacted powders from upper and lower layer of the sample. Whiskers
	from (a) the upper layer have relatively higher aspect ratio than from
	(b)the lower one.
	Synthesised powder from the lower layer of the sample	Synthesised powder from the upper layer of the sample
	SiC whiskers		Al 2 O 3 particles	SiC whiskers	Al 2 O 3 particles
	Length	Diameter	Size	Length	Diameter	Size
	10	0.5	6	20	0.6	4

FIG. 4. XRD patterns of heat treated kyanite-carbon powder at 1500 °C. The reaction for the formation of SiC and Al 2 O 3 is completed with the higher specific surface area carbon. (A) ␣-Al 2 O 3 ; (B) ␤-SiC; (M) mullite.