Maria Bekcheva

Cédric Join

Hugues Mounier

Hugues Mounier Cascaded

Cascaded model-free control for trajectory tracking of quadrotors

published or not. The documents may come

I. INTRODUCTION

In the last decade the use of unmanned aerial vehicles (UAVs) has increased significantly. From package delivery services, military uses and disaster management to photography and entertainment, the use-cases are numerous. Among the UAVs, quadrotors have been subject to intense research and development. One of the main challenges is modeling quadrotor dynamics and parameters estimation. This process can be tedious, time consuming and prone to modeling errors. A considerable amount of the literature has been devoted to model-based controls (Backstepping and sliding-mode techniques [START_REF] Bouabdallah | Backstepping and sliding-mode techniques applied to an indoor micro Quadrotor[END_REF], Inner-Outer Loop Control by applying PID controller for the inner-loop and nested saturation controller for the outer-loop [START_REF] Cao | Inner-Outer Loop Control for Quadrotor UAVs with Input and State Constraints[END_REF], Generalized Proportional Integral (GPI) observer based controller [START_REF] Sira-Ramirez | On the linear control of the quad-rotor system[END_REF], Model Predictive Control [START_REF] Alexis | Model predictive quadrotor control : attitude , altitude[END_REF]) among other methods. In these model-based approaches, the modeling accuracy can directly impact the controller success and performance. Moreover, many quadrotor applications are outdoor and may be faced with environmental uncertainties. Among the meteorological uncertainties, the wind effect has one of the highest impacts on the quadrotor performance while being very hard to predict. Therefore, in the model definition, the wind is usually considered as a predefined constant. As suggested in the papers [START_REF] Bouabdallah | PID vs LQ control techniques applied to an indoor micro quadrotor[END_REF], [START_REF] Pounds | Modelling and control of a large quadrotor robot[END_REF], [START_REF] Mellinger | Trajectory generation and control for precise aggressive maneuvers with quadrotors[END_REF], in practice, using a PID have shown more advantages than more advanced quadrotor controllers because of its simplicity and robustness. We chose to use the model-free control (MFC) since it has been already shown successful in practice (see [START_REF] Fliess | Model-free control[END_REF], [START_REF] D'andréa-Novel | A mathematical explanation via intelligent PID controllers of the strange ubiquity of PIDs[END_REF] and the references therein, and [START_REF] Abouaïssa | On ramp metering: towards a better understanding of ALINEA via model-free control[END_REF], [START_REF] Menhour | An Efficient Model-Free Setting for Longitudinal and Lateral Vehicle Control: Validation Through the Interconnected Pro-SiVIC/RTMaps Prototyping Platform[END_REF], [START_REF] Polack | Finite-Time Stabilization of Longitudinal Control for Autonomous Vehicles via a Model-Free Approach[END_REF]) that the MFC have better performances 1 Maria Bekcheva and Hugues Mounier are with the Laboratoire de Signaux et Systèmes (L2S, UMR CNRS 8506), CNRS, CentraleSupélec, Université Paris-Sud, France (e-mail: maria.bekcheva@l2s.centralesupelec.fr, hugues.mounier@l2s.centralesupelec.fr 2 Cédric Join is with CRAN (CNRS, UMR 7039), Université de Lorraine, BP 239, 54506 Vandoevre-lès-Nancy, France, A.L.I.E.N. (ALgèbre pour Identication & Estimation Numériques) and Projet Non-A, INRIA Lille -Nord-Europe, France (e-mail: cedric.join@univ-lorraine.fr) and advantages over a simple PID. These facts motivate us to present a cascaded-model-free approach for a quadrotor while considering more realistic situations like time-varying wind disturbances and measurement noise. We make our contribution in the following framework:

• We use a minimalist structure upon which our control scheme rests (see Equations ( 12), ( 14), (19a)-(19c)). This structure is independent of any mass, inertia, gyroscopic or aerodynamic effects; the only information we use is that each positional and rotational dynamics is of second order and that the thrust produced by the propellers is suitably projected onto the inertial frame. • Because of the underactuated nature of the quadrotor, we consider an inner-outer structure of the quadrotor that allows us to divide the control into two stages. At the first stage, we show that it is possible to apply a model-free control in the outer loop and from there, we derive the desired thrust, and the desired roll and pitch angles. We then apply again the model-free control in the inner loop for the second stage. • No precise information of the quadrotor physical parameters is required. For instance, it will be inconvenient and in some cases imprecise to measure the mass value [START_REF] Wang | Trajectory Tracking Control for Quadrotor Robot Subject to Payload Variation and Wind Gust Disturbance[END_REF] of the quadrotor together with its payload during each flight. • For the simulations, we have tested our approach with a realistic quadrotor model that contains gyroscopic and aerodynamic effects.

Recent studies have shown the interest in dividing the control problem into two parts: one based on the model and the other based on the model-free control to cope with model uncertainties and/or external perturbations. In [START_REF] Wang | Model-free-based terminal SMC of quadrotor attitude and position[END_REF], modelfree sliding-mode control (based on a linearized model) is studied, in [START_REF] Younes | Robust Model-Free Control Applied to a Quadrotor UAV[END_REF] a model-free backstepping control (based on a linearized model) and an LQR control have been experimentally tested on a real system. In these studies, the control law is still partially linked to the system model knowledge. In several works, the yaw angle is assumed to be zero due to high nonlinear coupling in position and attitude dynamics. In [START_REF] Wang | Event-driven model-free control in motion control with comparisons[END_REF], an event-driven model-free control is compared with several other controls while taking zero reference for the yaw angle. This limits the maneuverability of quadrotors. In our case we show that the quadrotor can follow aggressive maneuvers while tracking the yaw angle ψ. This is a huge advantage with clear benefits on time and energy. The positionyaw tracking as mentioned in [START_REF] Sanchez-Orta | Position-Yaw Tracking of Quadrotors[END_REF] can be useful in different use-cases e.g. Aerial Screwdriver task (the quadrotor should turn a screw), the Aerial Grasp (multiple quadrotors rotate their ψ-angle to grasp an object), etc. This paper has the following outline. In Section II, we briefly describe the quadrotor dynamics focusing on its structure. Then in Section III, we present our cascaded model-free approach. Finally, we validate our approach in Section V by testing it on two scenarios when the quadrotor is faced with unknown measurement noise and with unknown timevarying wind disturbances. The aggressive position trajectory is a B-spline curve with predefined control points.

II. QUADROTOR MODEL

The quadrotor is an underactuated system because it has six degrees of freedom but only four actual inputs. The six degrees of freedom include position motion in three directions and rotational motion around three axes. The schematic configuration of a quadrotor we adopted in this study is shown in Fig. 1 that includes the corresponding forces, angles and angular speeds.

A nonlinear model of quadrotor based on the Newton-Euler formalism (see the paper [START_REF] Luukkonen | Modelling and Ccontrol of Quadcopter[END_REF] for further information on the quadrotor model) is given by the following equations:

• Position dynamics mẍ = (sin ψ sin φ + cos ψ sin θ cos φ)u 1 -A x ẋ, (1a) mÿ = (-cos ψ sin φ + sin ψ sin θ cos φ)u 1 -A y ẏ, (1b) mz = -mg + cos θ cos φu 1 -A z ż, (1c)
• Attitude dynamics

I x φ = u 2 + θ ψ(I y -I z ) + J r θΩ r , (2a)
I y θ = u 3 + φ ψ(I z -I x ) -J r φΩ r , (2b)
I z ψ = u 4 + θ φ(I x -I y ), (2c)
where x, y and z are the position coordinates of the quadrotor's center of gravity, and θ, φ and ψ are the pitch, roll and yaw rotation angles respectively. The constant m is the mass, g is the gravitation acceleration, I x , I y , I z are the moments of inertia, and J r is the moment of inertia of the rotors. The controllers are: u 1 the total thrust generated by the four propellers applied in the z direction; u 2 , u 3 and u 4

Symbols and values Variables m = 0.53kg

Quadrotor mass

I x , I y = 6.228×10 -3 kg • m 2 Inertia parameters I z = 1.121×10 -2 kg • m 2 Inertia parameter J r = 6.01×10 -5 kg • m 2 Rotor inertia moment b = 3.13×10 -5 N • s 2 Thrust coefficient d = 7.5×10 -7 Nm • s 2 Drag coefficient A x , A y , A z = 0.25kg / s
Drag force coefficients the torques in the θ, φ and ψ directions respectively. The constants A x , A y and A z are the drag force coefficients for velocities in the corresponding directions of the inertial frame. The position dynamic model (1a)-(1c) contains the gravity force mg, the thrust force in the z direction and the drag force. The attitude dynamic model (2a)-( 2c) describing the roll, pitch and yaw rotations contains three terms which are the actuators action, the gyroscopic effect resulting from the rigid body rotation, and finally the gyroscopic effect resulting from the propeller rotation coupled with the body rotation. The control inputs u 1 , u 2 , u 3 , u 4 , and the speed Ω r are defined as:

u 1 = b(Ω 2 1 + Ω 2 2 + Ω 2 3 + Ω 2 4 ), u 2 = b(Ω 2 4 -Ω 2 2 ), u 3 = b(Ω 2 3 -Ω 2 1 ), u 4 = d(Ω 2 1 + Ω 2 3 -Ω 2 2 -Ω 2 4 ), Ω r = Ω 1 -Ω 2 + Ω 3 -Ω 4 (3)
where Ω i is the angular speed of the ith rotor, and b and d are the thrust and the drag coefficients respectively. This model will be used for the simulation. Table I lists the parameter values used for the simulation model in Section V.

III. CONTROL DESIGN

In this section, first we summarize some of the main ideas on Model-Free Control (MFC) introduced in [START_REF] Fliess | Model-free control[END_REF]. Then, we construct the cascaded-model-free controller for the quadrotor by using only a minimal nominal dynamics.

A. Preliminaries for Model-Free Control

The unknown differential equation describing the input/output behaviour of a finite-dimensional system with a single control variable u and a single output variable y E(y, ẏ, . . . , y (ι) , u, u, . . . , u (κ) ) = 0,

E : R ι+1 × R κ+1 → R, E ∈ C ∞ (R ι+1 × R κ+1 )
can be described as:

y (ν) = E(t, y, ẏ, . . . , y (ν-1) , y (ν+1) , . . . , y (ι) u, u, . . . , u (κ) ), (4)
where 0 < ν ι, ∂E ∂y (ν) = 0. For simplicity, the Equation (4) can be represented in short time interval by an ultra-local model as:

y (ν) = F + αu (5)
where

• α ∈ R is a non-physical constant parameter chosen such that αu and y (ν) will be of the same order of magnitude.

In the recent study [START_REF] Doublet | Model-free control for unknown delayed systems[END_REF], α is considered as a timevarying parameter in order to overcome the presence of unknown delays. • the time-varying function F (t) is approximated by a piecewise constant function. It represents the rest of the unmodelled dynamics in the input-output behaviour of the system and the unknown endogenous and exogenous disturbances. Therefore, it adapts to the changes of the system at each actuation step. Obtaining a good estimate of F can be achieved by considering it as ν th iterated integral on a short time interval [t-T, t] (see [START_REF] Fliess | Closed-loop parametric identification for continuous-time linear systems via new algebraic techniques[END_REF]). For instance, when ν = 2, we can rewrite the equivalent of the equation ( 5) in the Laplace form as:

s 2 Y (s) -sy(0) -ẏ(0) = F s + αU (s) (6)
We get rid of the initial conditions y(0) and ẏ(0) by differentiating twice the equation ( 6) w.r.t. to the Laplace variable s.

2Y (s) + 4s

dY ds + s 2 d 2 Y ds 2 = 2F s 3 + α d 2 U ds 2 (7)
To greatly attenuate the noise, we multiply the both sides of equation ( 7) by s -3

2 s 3 Y (s) + 4 1 s 2 dY ds + 1 s d 2 Y ds 2 = 2F s 6 + α s 3 d 2 U ds 2 (8)
In the time domain, we get

F = 5! 2T 5 t t-T (M 2 -4σM + σ 2 )y(σ) - α 2 M 2 σ 2 u(σ) dσ, (9)
where M = Tσ. The choice of the window T results in a trade-off. The larger is T , the smaller is the effect of noise but the larger is the error due to truncation. For a precise mathematical foundation on the treatment of the noise through iterated time integrals, see [START_REF] Fliess | Analyse non standard du bruit[END_REF], which is based on non-standard analysis, and also [START_REF] Mboup | Numerical differentiation with annihilators in noisy environment[END_REF] which rests on Jacobi polynomials. The closed-loop control applied to the ultra-local model (5) is defined as the so-called intelligent controller

u = - F -y (2)
d + C(e) α ( 10
)
where y d is the output reference trajectory, e = yy d is the tracking error, and C(e) = K p e + K d ė is a PD controller. Combining ( 5) and [START_REF] Fliess | Closed-loop parametric identification for continuous-time linear systems via new algebraic techniques[END_REF], it yields the following closed-loop error dynamics

ë + K d ė + K p e = e F = F -F . (11)
If the estimate F is good, the error e F 0 is small and choosing the gains such that K p > 0, K d > 0 guarantees a good tracking of y d .

B. Cascaded-model-free approach for the quadrotor Problem statement: We want to ensure that the quadrotor tracks the desired time-varying position trajectory (x d , y d , z d ) and yaw angle ψ d without precise information of the physical parameters and forces (for instance, the mass, the inertias and the aerodynamic forces) and despite external disturbances. We propose a cascaded model-free setting of the quadrotor based on a minimal nominal dynamics (see Equations ( 12), ( 14), (19a)-(19c)).

Inspired by the quadrotor dynamics structure, we divide the control in two stages (see Fig. 2): the outer loop with slow dynamics which controls the position, and the inner loop with fast dynamics which controls the attitude. From the quadrotor structure, we observe that the position dynamics depends on the attitude dynamics.

C. Outer-loop Position control

The model-free setting of the vertical z-dynamics (1c) is given by z = F z + α z u 1 [START_REF] Lyche | [END_REF] where F z represents the neglected dynamics and the external disturbances, and α z is a constant parameter. To estimate the time-varying F z over a time interval of length T , we use Equation [START_REF] Fliess | Model-free control[END_REF]. We set the control law for the thrust u 1 as:

u 1 = 1 α z (-Fz + zd -K pz e z -K dz ėz ) (13)
where e z = zz d and ėz = żżd are the tracking errors of the vertical position and velocity respectively. To operate safely the quadrotor requires a positive thrust u 1 > 0.

The model-free setting of the xy-dynamics (1a)-(1b) is given by

ẍ ÿ = F x F y + α xy R ψ sin φ sin θ cos φ u 1 . ( 14
)
where α xy = (α x , α y ) T are constants, and F x and F y are the neglected dynamics and disturbances for x and y respectively. The rotation matrix R ψ ∈ SO(2) is given by

R ψ = sin ψ cos ψ -cos ψ sin ψ , R -1 ψ = R T ψ , R ψ 1.
Here, • denotes the Euclidean norm.

Setting a new virtual input (v x , v y ) T defined as: we obtain

ẍ ÿ = v x v y = R ψ sin φ d sin θ d cos φ d u 1 , (15)
v x = 1 α x (-Fx + ẍd -K px e x -K dx ėx ), (16a)
v y = 1 α y (-Fy + ÿd -K py e y -K dy ėy ). ( 16b
)
where e x = xx d and e y = yy d are the tracking errors for x and y respectively. From there, we can deduce the reference trajectories θ d and φ d for the attitude dynamics by expressing the virtual input measures v x and v y in the body frame where ψ = 0 (no yaw dependence) by:

sin φ d sin θ d cos φ d = 1 u 1 R -1 ψ v x v y vx vy (17)
Hence we consider the following desired angles φ d and θ d :

φ d = arcsin(v x ), (18a)
θ d = arcsin vy cos φ d . (18b)

D. Inner-loop Attitude control

The attitude dynamics (2a)-(2c) are given by the following:

φ = F φ + α φ u 2 , (19a) θ = F θ + α θ u 3 , (19b) ψ = F ψ + α ψ u 4 (19c)
The attitude control performance is crucial since it is directly related to the actuators efficiency. After deducing the desired attitude φ d , θ d from the slow outer loop (18a)-(18b), the fast inner loop is stabilized by an MFC controller as

u 2 = 1 α φ (-Fφ + φd -K pφ e φ -K dφ ėφ ), (20)
u 3 = 1 α θ (-Fθ + θd -K pθ e θ -K dθ ėθ ), (21)
where e φ = φφ d , e θ = θθ d , ėφ = φφd , ėθ = θθd are the tracking errors of the angular positions and estimated angular velocities respectively. Similarly, for the yaw angle, we apply the Model-Free control as

u 4 = 1 α ψ (-Fψ + ψd -K pψ e ψ -K dψ ėψ ). (22)

IV. AGRESSIVE TRAJECTORY TRACKING

In this section, we give the predefined aggressive trajectory that the quadrotor will follow. In most cases, the trajectory is a straight line or circular orbit. We choose the B-spline curves (the definition of a B-Spline curve can be found in the Appendix) because they are always contained in the convex hull of their control polygon: the basis functions are positive and sum up to one (partition of unity), and have a local support [START_REF] Lyche | [END_REF]. By increasing the degree of the B-spline curve and/or by inserting extra knots, the distance between the control polygon and the B-spline curve can be reduced. We consider here an aggressive motion

[x d (t), y d (t), z d (t), ψ d (t)] : [t 0 t f ] → R 3 × SO(2)
defined as a B-spline curve:

x d (t) = 10 j=1 a xj B j,d , y d (t) = 10 j=1 a yj B j,d , z d (t) = 10 j=1 a zj B j,d , (23)
with degree d = 4, knot vector t = {0, 0, 0, 0, 0, 2.5, 5, 7.5, 10, 12.5, 15, 15, 15, 15, 15}, and control point vectors

A x = (a xj ) 10 xj=1 , A y = (a yj ) 10
yj=1
and A z = (a zj ) 10 zj=1 for x d (t), y d (t) and z d (t) respectively. We specialize the reference angle ψ d to a sigmoid between two constant angles a situation frequently needed in practice.

ψ d (t) = Ψ f -Ψ i 2 (1 + tanh(γ(t -t m ))) + Ψ i (24)
where Ψ i = 0rad is the initial ψ-angle and Ψ f = 1.5rad is the final ψ-angle of the quadrotor; γ = 0.4 is the slope parameter of the sigmoid and t m = 7 .

V. SIMULATION RESULTS

To demonstrate the proposed control we consider a realistic quadrotor simulator where we use the quadrotor model defined by the equations (1a)-(1c), (2a)-(2c) and (3) presented in Section II. Similar to the usual case in practice, the controller runs at 20Hz, and the sensor data is updated at 100Hz. The aim here is to track the positions x d , y d and z d , and the yaw angle ψ d . As initial positions and initial angles, we take (x 0 , y 0 , z 0 ) = (-1.8, 1, 0)[m] and (φ 0 , θ 0 , ψ 0 ) = (0, 0, 0)[rad] respectively. The initial linear velocities and initial angular velocities are zero. The control point vectors for the reference trajectories are the following:

P V P V P V K px 3 K dx 2 α x 1 K py 3 K dy 2 α y 1 K pz 10 K dz 15 α z 2 K pφ 2 K dφ 0.5 α φ 1 K pθ 1 K dθ 0.5 α θ 1 K pψ 3 K dψ 2.5 α ψ 9
A x = {-2, -1, 0, 2, 3, 3.5, 3.5, 5, 6.5, 7.5} ,
A y = {1.2, 2.5, 3.3, 1.8, 1.5, 2.5, 4, 4, 4, 4} , A z = {0, 0.2, 0.5, 2.5, 3.5, 3.9, 4.7, 6.5, 5.7, 3.5} .

(

) 25
Our trajectory is parametrized to avoid known static obstacles is a situation frequently needed in practice. In the nominal case (without disturbances), we observe in Fig. 3 (trajectory tracking), Fig. 4 (control inputs) and Fig. 5 (unknown dynamics estimation) satisfactory results when the position and the yaw angle change significantly. The 3D tracking is plotted in Fig. 6. The control gains that we consider are given in Table II. The values for α i are chosen from a large possible range each dependant on the corresponding state dynamics.

A. Scenario 1: Unknown measurement noise

For this scenario, we include the presence of the measurement noise b(t) ∼ N (0, σ 2 ) as an additive white Gaussian noise with zero mean and standard deviation σ = 0.15 in the four measured outputs x, y, z and ψ. The noisy measured outputs are set as:

x = x + b, ŷ = y + b, ẑ = z + b, ψ = ψ + b. (26)
Fig. 7 shows the noisy sensor measurements (green line) and the tracking trajectory of the quadrotor (blue line) that follows the desired trajectory (red line). From the result in Fig. 7 (trajectory tracking) and Fig. 8 (control inputs), it is apparent that the control approach is robust to the sensor disturbances without previous knowledge of the noise.

B. Scenario 2: Unknown time-varying wind disturbance

In this scenario, we investigate the quadrotor tracking in presence of wind disturbance that is not constant and is not assumed to be known. The high varying wind disturbances used in the simulations (displayed in Fig. 11) are represented by a sum of sinusoidal waves:

w(t) = 1.5µ 31 + µ 7 + 0.5µ 2 + 0.015µ 11 + 0.15b(t) (27)
where µ p (t) = sin tπ p and b(t) is the random Gaussian noise. We add this disturbance in the simulation model as an additive force disturbance such that mẍ = (sin ψ sin φ + cos ψ sin θ cos φ)u 1 -A x ẋ + w(t), (28a) mÿ = (-cos ψ sin φ + sin ψ sin θ cos φ)u 1 -A y ẏ + w(t), (28b) mz = -mg + (cos θ cos φ)u 1 -A z ż + w(t).

(28c)

We observe the unknowns F i for i = x, y, z, ψ in the nominal case (see Fig. 5) and in the presence of wind disturbance (see Fig. 10). We show that the precise estimation of the unknown F i for i = x, y, z, ψ leads us to a robust control performance. The trajectory tracked by the quadrotor in presence of time-varying wind disturbance depicted in Fig. 9 is almost identical to the nominal case.

VI. CONCLUSION

This methodology offers a cascaded model-free control design of a quadrotor. It stays robust, despite unknown disturbances. We have shown that a complex trajectory tracking on realistic scenarios is feasible. In the future works, we shall test our approach on an experimental platform. Considering the low complexity and the good performances of the proposed control, it sheds new lights on a possible commercial off-the-shelf solution. A complete stability analysis of the proposed system (quadrotor and cascaded modelfree controller) is in final stages and will be presented in a forthcoming publication. APPENDIX Definition 1: (B-splines [START_REF] Lyche | [END_REF]) Let the degree d be a nonnegative integer and let the knot vector t = (t j ), be a non-decreasing sequence of real numbers of at least d + 2. The jth B-spline of degree d with knots t is defined by:

B j,k,t (x) = x -t j t j+d -t j B j,d-1,t + t j+k -x t j+1+d -t j+1 B j+1,d-1,t (29)
for all real numbers x, with

B j,0,t (x) = 1, if t j x < t j+1 0, otherwise. (30)
Definition 2: (B-spline curves). Let a knot vector for a total of n B-splines be defined as:

t = (t j ) m=n+d+1 j=1 = [0, . . . , 0 d + 1 , t d+1 , . . . , t m-d-1 , 1, . . . , 1 d + 1 ]
The linear space of all linear combinations of these B-splines is the spline space S d,t defined by

S d,t = n j=1 c j B j,d | c j ∈ R for 1 j n

Fig. 1 .

 1 Fig. 1. The quadrotor system.

u 1 θ d ; φ d u 2 Fig. 2 .

 122 Fig. 2. Cascaded control overview (The subscript d denotes desired trajectory.)

Fig. 3 .

 3 Fig.3. Reference trajectories and model-free tracking trajectories of the four outputs x, y, z and ψ in the nominal case

Fig. 4 .Fig. 5 .Fig. 7 .Fig. 8 .

 4578 Fig. 4. Control inputs in the nominal case

Fig. 9 .Fig. 10 .Fig. 11 .

 91011 Fig. 9. Tracking outputs in presence of wind disturbance

TABLE I THE

 I PARAMETERS AND THEIR CORRESPONDING VALUES FOR THE QUADROTOR SIMULATION.


Adnane Hamiaz
email: ahamiaz@taibahu.edu.sa

Xavier Ferrieres
email: xavier.ferrieres@onera.fr

Olivier Pascal
email: olivier.pascal@laplace.univ-tlse.fr

Efficient Numerical Algorithm to Simulate a 3D Coupled Maxwell-Plasma Problem

Keywords: Maxwell's equations, fluid plasma model, time-domain

This paper proposes an improved algorithm, based upon an explicit finite difference scheme, in order to simulate the plasma breakdown induced by a monochromatic High Power Micro-Wave (HPM). The 3D coupled Maxwell-plasma equations are to be solved. We want to study with this model the geometry of the discharge and plasma formation at high pressure which may contribute to shield microwave sensors or circuits. Generally, the simulation of this kind of problem is very time-consuming, but by using the fact that the plasma evolution in time is slow relatively to the monochromatic source period, we can drastically reduce the simulation time. By considering this assumption, we describe in the paper a process which allows to obtain this important reduction. Finally, an example where we show the gain obtained in terms of computation time with our process is given to validate and illustrate the global work.

simulations, Yee's scheme and finite-difference methods.

Introduction

In this paper, we are interested in the breakdown induced by an electromagnetic HPM (High Power Micro-Wave). An efficient numerical method is proposed, based upon a finite differences scheme, to reduce the computation time in the simulation of a 3D coupled Maxwell-Plasma problem. In particular, we are interested in studying the microwave breakdowns on small structures and plasma formation at high pressure.

Such a kind of problem presents a big interest in ElectroMagnetic Compatibility (EMC) for the evaluation of the electromagnetic fields induced by an external source, on electronic components located inside a cavity which is not closed. By considering a HPM source, the important value of the electric field can create discharges. They are located on the small apertures or on the thin slots of the cavity. They are able to protect the inside components by reducing the value of the electromagnetic fields. The absorbing and reflecting properties of the plasma medium then enhance the shielding but it can also deteriorate dielectric materials in the long time as it has been observed in the context of partial discharges.

In the first part of the paper, we recall the equations of the studied physical problem proposed by J.P.Boeuf to solve the plasma formation at high pressure in a 2D context [START_REF] Boeuf | Theory and modelling of selforganization and propagation of filamentary plasma arrays in microwave breakdown at atmospheric pressure[END_REF]. In the second part, a 3D approximation to simulate this physical problem is detailed. This approximation, already developed in a 2D context [START_REF] Boeuf | Theory and modelling of selforganization and propagation of filamentary plasma arrays in microwave breakdown at atmospheric pressure[END_REF] [START_REF] Chaudhury | Pattern formation and propagation during microwave breakdown[END_REF], is based on an explicit finite differences scheme and presents an important computational cost due to the different time and space characteristics of the electromagnetic problem and of the plasma problem. Indeed, to ensure an accurate solution for the plasma problem, it is mandatory to use a very small step size in space which implies for the Maxwell problem a very small time step to ensure stability. In the 3D case, this drawback involves an excessive computing time to treat efficiently the entire problem. To avoid this drawback, some authors have proposed an implicit numerical scheme to solve the Maxwell problem [3] [4], but this solution don't take into account the great diversity of physical model available in the Yee scheme for thin wires, slots, ..., needed in many of our applications. Therefore, we propose an optimized explicit numerical scheme where the new idea is to reduce the cost of the electromagnetic fields computation by avoiding to evaluate these fields at each step of the scheme. A quasi stationary state is assumed for the plasma medium to extract the fields in an appropriate number of periods. This allows to update the plasma density with a much longer time step without calculating electromagnetic fields all along.

In the third part of this paper, we give the principle of this optimization and an example to validate this approach and evaluate the gain in terms of computation time.

Physical model

The mathematical formulation of the physical model, taken into account for our study, is based upon a coupling between Maxwell's equations and plasma's equations.

In the literature, a lot of papers are dealing about the resolution of this system of equations. In particular, we can cite, the works of Yu [START_REF] Yu | A 3 3-d global earth-ionosphere fdtd model including an anisotropic magnetized plasma ionosphere[END_REF] and Simpson [START_REF] Samini | An efficient 3-d fdtd model of electromagnetic wave propagation in magnetized plasma[END_REF] for studying the electromagnetic waves propagation inside a magnetized plasma. However, concerning our applications with a non magnetized medium, the most interesting model is the one of Boeuf [START_REF] Boeuf | Theory and modelling of selforganization and propagation of filamentary plasma arrays in microwave breakdown at atmospheric pressure[END_REF][2] [START_REF] Chaudhury | Computational study of filamentary pattern formation in a high power microwave breakdown generated air plasma[END_REF], where the mathematical formulation is given by :

µ 0 ∂H ∂t + ∇ × E = 0 0 ∂E ∂t -∇ × H = -J J quv (1)
where (E, H) are the electromagnetic fields, J the current density, and e, q = -e, u, v are respectively the elementary charge, the electron charge, the density and the mean velocity. The electron velocity is obtained from a simplified momentum transfer equation:

∂v ∂t + ν c v = q m e E (2)
where m e is the electron mass, ν c = 5.3 × 10 9 p r is the electron-neutral collision frequency in air and p r is the pressure in torr.

The plasma is assumed to be quasi-neutral (n e = n i = u) and the time evolution of the plasma density u is described by the diffusion equation :

∂u ∂t -∇ • (D ef f ∇u) = ν ef f u (3)
where ν ef f = ν i -ν a is the effective ionization frequency taking into account ionization (ν i ) and attachment (ν a ) frequency. It was shown in [START_REF] Chaudhury | Pattern formation and propagation during microwave breakdown[END_REF] that the more exact conservative form of the diffusion term in equation ( 3) leads to very similar numerical results. The values for ν i and ν a are functions of local effective field E ef f and are obtained in air from BOLSIG+ [START_REF] Hagelaar | Solving the boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models[END_REF]. We can also find the expression of ν ef f in the report [START_REF] Larry K Warme | Ionization coefficient approach to modeling breakdown in nonuniform geometries[END_REF]. The local effective field can be written as :

E ef f = E rms 1 + ω 2 /ν 2 c ( 4
)
where E rms is the local root mean square of the electrical field, and ω is the angular frequency of the incident field.

D ef f is an effective diffusion coefficient defined as

D ef f αD e + D a 1 + α , α = ν i ε 0 q u µ e
where ν i is the ionization frequency, D e = µ e k B T e /q is the free electron diffusion coefficient and D a = µ i k B T e /q is the ambipolar diffusion coefficient. T e is the electron temperature (in Kelvin), µ e = q/(m e ν c ) is the electron mobility, µ i is the ion mobility and k B is the Boltzmann constant. It was shown heuristically in [START_REF] Boeuf | Theory and modelling of selforganization and propagation of filamentary plasma arrays in microwave breakdown at atmospheric pressure[END_REF] and confirmed numerically in [START_REF] Chaudhury | Computational study of filamentary pattern formation in a high power microwave breakdown generated air plasma[END_REF] that the effective diffusion coefficient above can satisfactorily describe the fact that the plasma front, where the electron and ion densities are small and the Debye length is large, diffuses with the free electron diffusion coefficient while the plasma bulk diffusion is ambipolar. We assume in the simulations that k B T e is constant and equal to 2 eV , and that µ i = 10 -2 µ e (so that D a = 10 -2 D e ) [START_REF] Chaudhury | Pattern formation and propagation during microwave breakdown[END_REF].

For more details about this physical model, we can see the PhD thesis of Guo Quang Zhu [START_REF] Zhu | Modeling of plasma dynamics and patterne forulation during high pressure microwave breakdown[END_REF] and the PHD thesis of K. Kourtzanidis [START_REF] Koutzanidis | Numerical simulation of plasma acuators for flow control[END_REF].

FDTD approximation

To solve Maxwell's equations in time domain, a finite-difference method (FDTD) was developed in 1966 by K. S. Yee [START_REF] Yee | Numerical solution of initial boundary values problems involving maxwells equations in isotropic media[END_REF]. It is currently one of the most famous method used to solve these equations. To solve our coupled Maxwellplasma problem, we are going to adapt the Yee scheme to our system of equations formed by ( 1), ( 2), (3) :

                           0 ∂E ∂t -∇ × H = -J µ 0 ∂H ∂t + ∇ × E = 0 J = quv ∂u ∂t -div(D ef f ∇u) = ν ef f u ∂v ∂t + ν c v = q me E u = 0 (x, y) ∈ ∂Ω (5)
The Yee's scheme is based on a leap-frog scheme in both space and time. We define the computation domain Ω as a parallelepiped where some absorbing boundary conditions are applied to emulate the infinite space [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF].

In our case, we consider an uniform cartesian mesh and we split Ω such that :

Ω = [0, Lx] × [0, Ly] × [0, Lz] = i,j,k Ω i,j,k
where

Ω i,j,k = [x i , x i+1 ] × [y j , y j+1 ] × [z k , z k+1 ] so that x i = (i -1)∆x, y j = (j -1)∆y, z k = (k -1)∆z, for 1 ≤ i ≤ N x , 1 ≤ j ≤ N y , 1 ≤ k ≤ N z , where ∆x = Lx/N x , ∆y = Ly/N y , ∆z = Lz/N z and N x , N y , N z ∈ N * .

Approximation of Maxwell's Equations

75

The discretization of the Maxwell equations by the Yee scheme is well known and has been developed in many papers. For a complete study of this scheme, we can see, for example, the book of Taflove [START_REF] Taflove | Computational Electrodynamics: The Finite-Difference Time-Domain Method[END_REF]. Here and for an easier understanding, we recall that the Yee scheme is used to simulate the following system of equations :

   ε ∂E ∂t -∇ × H = -J µ 0 ∂H ∂t + ∇ × E = 0 (6)
where the unknowns (E x , E y , E z ) and (H x , H y , H z ) of the scheme are located as shown in Figure ( 1) and evaluated at time t = n∆t for the electric fields and t = (n + 1/2)∆t for the magnetic fields. More precisely, the components evaluated by the scheme are given by :

E n i,j,k =       E n x i+ 1 2 ,j,k E n y i,j+ 1 2 ,k E n z i,j,k+ 1 2       , H n+1/2 i,j,k =      H n+1/2 x i,j+ 1 2 ,k+ 1 2 H n+1/2 y i+ 1 2 ,j,k+ 1 2 H n+1/2 z i+ 1 2 ,j+ 1 2 ,k     
where, for a given component U , U n i,j,k defines the value of this component taken at the time t = n∆t and at the location x = i∆x, y = j∆y, z = k∆z, where ∆t, ∆x, ∆y, ∆z define respectively the steps in time and space domains. Finally, the Yee scheme is explicit and of two order in time and in space when homogeneous step sizes are considered. The well-known drawbacks of this scheme are the staircase representation of the objects and its significant numerical dispersive error. In free space, the scheme is stable under the condition [START_REF] Taflove | Computational Electrodynamics: The Finite-Difference Time-Domain Method[END_REF] :

∆t ≤ 1 ν 1 1 ∆x 2 + 1 ∆y 2 + 1 ∆z 2 (7)
where ν = 1/ √ ε 0 µ 0 defines the velocity of the waves inside the medium.

Approximation of the plasma equations 80

Consider the plasma's equations given in the physical model by :

   ∂u ∂t -∇ • (D ef f ∇u) = ν ef f u ∂v ∂t + ν c v = q me E (8)
In our finite differences approximation, the unknowns are the plasma density u and the velocity v = (v x , v y , v z ). By considering the same mesh as for the electromagnetic fields, we choose to locate the plasma density at the mesh nodes while the velocity position is identical to the one of electric field in the Yee scheme (see Figure 2). By using the same notation for the plasma unknowns as for the electromagnetic fields, the global finite differences discretization for the velocity equation is given by :

v n+1/2 i,j,k -v n-1/2 i,j,k ∆t + ν c v n+1/2 i,j,k + v n-1/2 i,j,k 2 = q m e E n i,j,k (9)
with

v n+1/2 i,j,k =      v n+1/2 x i+ 1 2 ,j,k v n+1/2 y i,j+ 1 2 ,k v n+1/2 z i,j,k+ 1 2     
and for the diffusion equation which describes the plasma density, we have :

u n+1 i,j,k -u n i,j,k ∆t -(∇ • (D ef f • ∇u)) n i,j,k = ν ef f i,j,k u n i,j,k (10)
We reformulate the operator (∇ • (D ef f ∇u)) as the sum of two operators:

∇ • (D ef f ∇u) = ∇D ef f • ∇u + D ef f ∆u
The first term on the right side of the previous equation is discretized as follows :

(∇D e • ∇u) n i,j,k =

D n ef f i+1,j,k -D n ef f i,j,k x i+1 -x i u n i+1,j,k -u n i,j,k x i+1 -x i + D n ef f i,j+1,k -D n ef f i,j,k y j+1 -y j u n i,j+1,k -u n i,j,k y j+1 -y j + D n ef f i,j,k+1 -D n ef f i,j,k z k+1 -z k u n i,j,k+1 -u n i,j,k z k+1 -z k (11)
and the discretization of the second term is given by :

(D ef f ∆u) n i,j,k = D n ef f i,j,k [ 2(u n i+1,j,k -u n i,j,k ) (x i+1 -x i )(x i+1 -x i-1 ) + 2(u n i-1,j,k -u n i,j,k ) (x i -x i-1 )(x i+1 -x i-1 ) + 2(u n i,j+1,k -u n i,j,k ) (y j+1 -y j )(y j+1 -y j-1 ) + 2(u n i,j-1,k -u n i,j,k ) (y j -y j-1 )(y j+1 -y j-1 ) + 2(u n i,j,k+1 -u n i,j,k ) (z k+1 -z k )(z k+1 -z k-1 ) + 2(u n i,j,k-1 -u n i,j,k ) (z k -z k-1 )(z k+1 -z k-1 ) ] (12)
For the effective diffusion coefficient D ef f , we take :

D n ef f i,j,k = D a + α n i,j,k D e 1 + α n i,j,k with α n i,j,k = ν i ε 0 /(q u n i,j,k µ e ) (13)
In the case where D ef f is constant and ν ef f = 0, the scheme ( 10) is stable under the condition (see Appendix A) :

∆t ≤ 1 2D e 1 1 ∆x 2 + 1 ∆y 2 + 1 ∆z 2 (14)
with (D a ≤ D ef f ≤ D e ).

Time discretization of the global problem and stability :

In the numerical solutions of equations ( 5) two different time steps ∆t M and ∆t u are used respectively for the discretization of the Maxwell and plasma equations. This is due to the fact that the time evolution of the plasma density is very slow compared to that of the electromagnetic field which is related to the incident wave frequency.

In the equations ( 5), we choose to take an explicit Euler scheme, a leapfrog scheme and an implicit centred scheme to discretize in time respectively the plasma density equation, the Maxwell equations and the velocity equation.

In our time approximations, we consider the second members in both Maxwell and plasma equations (5) taken at the previous time. By considering this assumption, the stability studies are independently done on every equation of the system (5) and we find two conditions on the time steps :

     ∆t M ≤ 1 ν 1 1 ∆x 2 + 1 ∆y 2 + 1 ∆z 2 ∆t u ≤ 1 2De 1 1 ∆x 2 + 1 ∆y 2 + 1 ∆z 2 (15)
By considering our numerical experiments, theses conditions seem sufficient to ensure the stability of the numerical scheme. The global theoretical condition for the stability of our coupled systems of equations is difficult and remains an open mathematical problem which is not treated here. In our simulations, the values of D e , ν and of the minimum cell dimension are such that the ratio between ∆t u and ∆t M is quite large. Indeed, by considering an homogeneous mesh where ∆ x = ∆ y = ∆ z = dl, we have ∆t M ∆tu = De 2νdl . In this formula, the velocity of the waves in the medium is equal to ν = 3.e8ms, the coefficient D e , after calculation by considering the values previously given for k B T e /q, nu c and nu e = q/(m e ν e ), is lower than 1 and the step taken into account in our simulations dl = 1e -4m.

With these values, the ratio ∆t M ∆tu is of the order of 3000. Then to reduce the computation time in our scheme, we use two time steps as in [START_REF] Chaudhury | Computational study of filamentary pattern formation in a high power microwave breakdown generated air plasma[END_REF][15] [START_REF] Samini | An efficient 3-d fdtd model of electromagnetic wave propagation in magnetized plasma[END_REF]. ∆t M to compute the electromagnetic fields and ∆t u = m ∆t M to evaluate plasma density, where m is an integer value given by m = int(

∆t u ∆t M )
1, where ∆t u is obtained by the stability condition for the density plasma. In the calculation process of the electromagnetic fields, between the two steps ∆t M and m ∆t M , the plasma density is assumed constant and equal to the value evaluated at ∆t M .

Concerning the velocity equation, the time step taken into account is the same as for the electromagnetic fields.

Finally the numerical scheme for the system (5), is given by :

                     E n+1 i,j,k = E n i,j,k + ∆t M ε0 ∇ h × H n+1/2 -J n+1/2 i,j,k H n+1/2 i,j,k = H n-1/2 i,j,k -∆t M µ0 ∇ e × E n v n+1/2 i,j,k = (1 -∆t M νc 2 )v n-1/2 i,j,k + ∆t M q me E n i,j,k /(1 + ∆t M νc 2 ) u N +1 i,j,k = u N i,j,k + ∆t u [(∇ • (D ef f • ∇u)) N i,j,k + ν ef f i,j,k u N i,j,k ] J n+1/2 i,j,k = q u N -1 i,j,k +u N i,j,k 2 v n+1/2 i,j,k (16)
with t n = n∆t M and t N = N ∆t u = N m ∆t M ,

Optimization of the simulation time

The use of two different time steps to evaluate the electromagnetic fields and 95 the plasma density allows to reduce the computation time, however it remains important. This is the reason why we propose a new idea to reduce the number of evaluations of the electromagnetic field by considering the monochromatic property of the source and the relatively slow plasma density evolution.

In this section, we present this new algorithm based upon the facts that : • the plasma density is constant between t N and t N +1 and equal to u N ,

• the plasma density depends only of the E rms value and is then not instantaneously dependent of the electromagnetic field (E, H),

• the electric and magnetic fields are periodic in time due to the monochromatic source.

Considering an incident monochromatic field with a period T and assuming that on the computational domain Ω the electromagnetic fields have a stable permanent regime after p periods (we call p the parameter of optimisation).

Let t N be a given time where we have computed a plasma density u N and an electromagnetic field (E N , H N ). We want to compute these quantities at the time t N +1 = t N + ∆t u , where ∆t u is the time step to evaluate the plasma density. In order to optimize the computation time, the idea is to perform the following steps (see Fig. 3).

• First, we re-evaluate ∆t u so that ∆t u = kpT , where k is defined by the floor function value of ∆t u /(pT ) ;

• Evaluate the electromagnetic fields on the interval [t N , t N + pT [, with t N + pT t N +1 , on the different time points given by the time step ∆t M , related to the Maxwell equations. In this step we consider the plasma density constant and equal to u N ;

• Evaluate the RMS field E rms (x, y, z) = 1 T t N +pT t N +(p-1)T E(t, x, y, z) 2 dt ;

• Evaluate the plasma density at the time step t N +1 ;

• Evaluate the electromagnetic fields at t N +1 . By using the periodicity, the electromagnetic fields values at t N +1 and at t N + pT are the same. Then, there are no calculation to do for evaluating the electromagnetic fields at implies an important reduction of the computation time, if ∆t u is very larger than pT . This last condition is generally true in our applications.

Numerical results

In this section, we present some results on the gain, in terms of computation time, obtained by our optimized algorithm for simulations on the Maxwell problem and the Maxwell-plasma problem. Concerning the Maxwell-plasma problem, we propose to study the evolution of a plasma density of a discharge at high pressure in air. To make this simulation, we consider a plane wave which illuminates an initial ball of low density plasma, in order to localize the plasma discharge in space without any focusing. The initial density is low enough that it remains transparent to the wave energy but able to later generate the breakdown.

Description of the configuration

The configuration studied consists in the scattering of a plane wave on a plasma discharge in a 3D infinite domain. We evaluate for this configuration, Maxwell problem and the following parameters to solve the plasma equations :

• a time step for the plasma ∆t u = 1.e -9s = 100T , where T definies the period of the incident plane wave signal ;

• ν m = 3.8e12s -1 .

180

• an initial plasma density (m -3 )

u 0 (x, y, z) = 1.e18 exp( -((x -x 0 ) 2 + (y -y 0 ) 2 + (z -z 0 ) 2 ) 2σ 2
)), with σ = 1.e -5m and P 0 = (x 0 , y 0 , z 0 ) = (4.5mm, 1.5mm, 1.5mm).

The Figure [START_REF] Yu | A 3 3-d global earth-ionosphere fdtd model including an anisotropic magnetized plasma ionosphere[END_REF] shows the evolution of the electric field at the point P 1 obtained with and without the optimized algorithm. We can again observe a good agree-ment between both results. Concerning the plasma density, Figure [START_REF] Samini | An efficient 3-d fdtd model of electromagnetic wave propagation in magnetized plasma[END_REF] shows the value obtained at the test-point P 1 for both approaches and the differences between these solutions obtained with and without our optimized algorithm.

We note on this figure, the very good agreement between both solutions with a same gain in terms of computation time regarding the Maxwell simulation, for our optimized algorithm.

5.4. Improving the parameter of optimization p in the evaluation of the electromagnetic fields

In the previous example, using the optimized algorithm, we evaluated the electromagnetic fields during a time interval equal to ten periods (p = 10) of the incident plane wave. In this subsection, we show that it is possible to decrease the computation time by reducing this time interval to three periods (p = 3), without modifying the accuracy of the solution. By taking into account the previous example where we carried out a Maxwell-plasma simulation, we make different computations by choosing a parameter of optimization p equal to 3, 4 and 5. This parameter is used to calculate the time interval in which we solve Maxwell's equations and the velocity's equations. Figure 7 and table 1 present the plasma density at the test-point P 1 and the computation time consumed for different values of p. By considering these results we notice that we can improve the computation time given in the previous example by a factor 3. In our configuration, it was not possible to take a time interval less than 3 periods to evaluate the electromagnetic fields because they are not in a stable state. This parameter appears to be physically dependent on the non resonant behaviour of the presented configuration and the travelling time of the incident wave through the domain. It is then clearly related to each physical problem under investigation. The validation cases proposed here are therefore very favorable regarding the benefit of the optimized process. Indeed, for the study of discharges and plasma formation at atmospheric pressure due to High Power Microwave, in our applications the field amplitude and the source frequency considered are respectively greater than 3M V and 1GHz. With these assumptions, before having a modification of the density plasma by diffusion, the electromagnetic field travels a significant number p of periods (see [START_REF] Guo Qiang Zhu | Ionization-diffusion plasma front propagation in a microwave field[END_REF]). To validate and to quantify the performances of our approach on a physical problem describing the plasma formation at atmospheric pressure, we propose to compared it with some results obtained in [START_REF] Kourtzanidis | Three dimensional simulations of pattern formation during high-pressure, freely localized microwave breakdown in air[END_REF] with an implicit scheme to solve Maxwell-plasma equations. To have similar simulations in terms of memory and computational costs, we modify our previous configuration, by taking a source frequency equal to f 0 = 110GHz, and a domain of interest defined by 2λ × 2.2λ × 1.5λ, where λ is the wavelength of the source. With this new configuration, we make simulations by using 64 processors and an approximation of the computational domain given by 400×440×300 cells, as for [START_REF] Kourtzanidis | Three dimensional simulations of pattern formation during high-pressure, freely localized microwave breakdown in air[END_REF]. Figure [START_REF] Hagelaar | Solving the boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models[END_REF] 225

shows the evolution of the density plasma at different time-step, obtained with our approach. We note on this figure, the same behaviour of the density plasma as in the results given in the paper [START_REF] Kourtzanidis | Three dimensional simulations of pattern formation during high-pressure, freely localized microwave breakdown in air[END_REF], with moreover, an important gain in computation time. Indeed, the time-consuming with our optimized explicit method approach is approximatively equal to ten hours, against one week for 230 the implicit scheme.

Conclusion

In this paper, we proposed and validated an optimized algorithm to efficiently compute electromagnetic fields and plasma densities in a coupled Maxwell-Plasma problem where the source is a monochromatic HPM that ignites the breakdown. This new algorithm permits to obtain for 3D simulations an important gain in terms of computation time. This gain is higher when dealing with limited size domain, in terms of wavelengths, and non, or low, resonant cases. Because using the conventional FDTD method without this optimization is difficult and sometimes impossible to properly handle a 3D Maxwell-plasma problem, the proposed process is of great interest. With our optimized algorithm, it is now possible to hope handling larger problems and investigating specific 3D non linear evolutions of plasmas and electromagnetic fields.

with A = cos(kx∆x)-

  70

Figure 1 :

 1 Figure 1: Localisation of the unknown in then FDTD approximation.

  85

Figure 2 :

 2 Figure 2: Unknown position related to the plasma in the cell grid

100•

  we have two different time steps ∆t M and ∆t u for the discretization of the Maxwell and plasma equations,

Figure 3 :

 3 Figure 3: Optimized algorithm.

Figure 4 :

 4 Figure 4: Evolution of the Ey component of the electric field (top: left) at the point P 1 with and without optimized algorithm. Other figures are the zoom at different times.

Figure 5 :

 5 Figure 5: Evolution of the Ey component of the electric field at the test-point P 1 with and without optimized algorithm for different ranges of time.

  215

Figure 6 :

 6 Figure 6: Evolution of the plasma density (left) and the relative error (right) at the test-point P 1 with and without optimized algorithm.

Figure 7 :

 7 Figure 7: Evolution of the plasma density at the test-point P 1 with and without optimized algorithm by changing the parameter p.

Figure 8 :

 8 Figure 8: 3D-spatial evolution of the electron density isosurface (1.e21m -3 ) contours coloured by the electric fields, at 19ns, 27ns, 38ns, 50ns, 59ns and 70ns.

Table 1 :

 1 Computation time as a function of p 5.5. Evolution of the spatial density plasma in a discharge in air at high pressure

Acknowledgement

The first author thanks Taibah university for her research funding. This is an extension of a previous work funded by PLASMAX project.

the spatial evolution of the plasma density and the electromagnetic fields along the time.

First, to study the performance of our optimized method, we consider a rectangular computational domain defined by [0, 6mm] × [0, 3mm] × [0, 3mm] and bounded by Perfectly Matched Layers (PML) to simulate the infinite space.

A monochromatic (f 0 = 100GHz) incident plane wave source is defined by

where at a point (x, y, z) and at the time t the electric field is given by E y (x, y, z) = E 0 sin(ω(t -k x • x/c 0 )), with E 0 = 6M V /m, c 0 is the speed of the wave in the vacuum and ω = 2πf 0 . The initial plasma density is localized inside a ball of center P 0 = (4.5mm, 1.5mm, 1.5mm) and with a radius r = 0.1mm. The initial value of the plasma density is equal to u 0 = 1.e18m -3 .

The problem is considered at atmospheric pressure (pr = 760T orr) and we want to follow the evolution of the electromagnetic fields and the plasma density at the test-point P 1 = (4.5mm, 1.5mm, 1.5mm). In our first simulations, we have fixed the optimisation parameter p introduced in the previous section to 10.

Next, we propose to evaluate with our approach, 3D spatial plasma density evolution on an example proposed in [START_REF] Kourtzanidis | Three dimensional simulations of pattern formation during high-pressure, freely localized microwave breakdown in air[END_REF] and compare the computation time for the two methods.

Analysis of performance : Maxwell simulations

Considering the previous configuration, first, we show the advantage of using our optimized algorithm to evaluate only the electromagnetic fields by assuming no plasma density in the problem. In our example, we take a time step ∆t M = 2.5e -14s. Figure [START_REF] Berenger | An implicit fdtd scheme for the propagation of vlf-lf radio waves in the earth-ionosphere waveguide[END_REF] shows the fields computed by using our optimized algorithm and the non optimized method. We can see that the fields computed with both methods are in very good agreement. The great interest of our optimized algorithm is the reduction of 90% of the computation time to have a similar solution than the non optimized method.

Analysis of performance : Maxwell-plasma simulation

Now, in this section, we are interested in the Maxwell-plasma problem. In our simulations, we take the same time step ∆t M as previously to solve the Appendix A By considering D ef f = D e constant and ν ef f = 0, the equation ( 10) is given by :

with u n i,j,k = U (t = n∆t, x i = i∆x, y j = j∆y, z k = k∆x). To establish a condition of stability, we use the Von Neumann stability analysis by introducing in the equation ( 17) a modal solution given by u n i,j,k = u n e i(kxxi+kyyj +kzz k ) . We obtain :

which can be written :