E Nardon

A Fil

P Tamain

R Guirlet

H R Koslowski

M Lehnen

C Reux

F Saint-Laurent

On the mechanisms governing gas penetration into a tokamak plasma during a massive gas injection

A new 1D radial fluid code, IMAGINE, is used to simulate the penetration of gas into a tokamak plasma during a Massive Gas Injection (MGI). The main result is that the gas is in general strongly braked as it reaches the plasma, due to mechanisms related to charge exchange and (to a smaller extent) recombination. As a result, only a fraction of the gas penetrates into the plasma. Also, a shock wave is created in the gas which propagates away from the plasma, braking and compressing the incoming gas. Simulation results are quantitatively consistent, at least in terms of orders of magnitude, with experimental data for a D 2 MGI into a JET Ohmic plasma. Simulations of MGI into the background plasma surrounding a runaway electrons beam show that if the background electron density is too high, the gas may not penetrate, offering a possible explanation to the recent results of Reux et al. in JET [1].

I. INTRODUCTION

Due to the potentially deleterious consequences of disruptions in ITER, a Disruption Mitigation System (DMS) is needed [START_REF]on Disruptions, Plasma Control, and MHD and ITER Physics Basis Editors[END_REF] [3]. Its design, which is presently underway, is based on the use of massive material injection, in the form of either gas (Massive Gas Injection or MGI) or shattered pellets [4] [5]. This paper focuses on MGI and more specifically on its modelling. Our objective is to identify key mechanisms which determine the penetration of the gas into the plasma and to study their dependency on plasma parameters. MGI consists in releasing a large amount of gas (compared to the plasma content) by suddenly opening a gas reservoir which can be located from a few cm to a few m away from the plasma. After a certain time of propagation in vacuum, the arrival of the gas into the plasma typically leads to an increase of the electron density n e and the propagation of a cold and radiating front into the plasma (see for example Figure 12 in [6] for high-resolution profiles in MAST). At some point, a thermal quench (TQ) is triggered. There is evidence in several machines (e.g. MAST [6], Tore Supra [8] and TEXTOR [7]) that the TQ is triggered when the cold front reaches the q=2 surface.

Significant work has been devoted by experimentalists to the estimation of the ratio between the increase in the number of electrons ∆N e and the number of gas atoms that has reached the edge of the plasma, N inj (the latter being estimated using models for the gas propagation into vacuum) [8] [9] [10] [11]. This ratio, which is often called "fuelling efficiency" and noted F ef f , is indeed of importance regarding, for example, the objective of suppressing Runaway Electrons (RE) by increasing n e [4] [5]. Depending on the type of gas, the reservoir pressure and characteristics (e.g. its distance to the plasma) and the plasma parameters, F ef f can vary significantly, typically from a few % to a few tens of %. This indicates that the gas penetrates only partially, and sometimes even marginally, into the plasma. This is consistent with, for example, the observation from DIII-D that "over a wide range of initial target conditions [...], visible camera images of neutral Ar emission indicate that the propagation of jet neutrals is stopped at the plasma edge ( 0 -5 cm past the separatrix) during the TQ" [START_REF] Hollmann | DIII-D Studies of Massive Gas Injection Fast Shutdowns for Disruption Mitigation[END_REF]. Another important recent result possibly related to the question of gas penetration is the unsuccessful attempt to suppress an RE beam after its formation using MGI on JET [1]. One possible reason for this failure is indeed that the gas does not penetrate deep enough into the background plasma and thus cannot reach the RE beam. Clearly, the above points motivate efforts to understand the mechanisms governing gas penetration during an MGI. The development of a new 1D radial fluid code called IMAGINE and its application to MGI simulations described below are part of these efforts.

A number of codes have been used for MGI modelling over the past, including ASTRA [START_REF] Leonov | [END_REF] [14], TOKES [15] [16], SOLPS [START_REF] Pautasso | Proc. 22nd Int. Conf. on Fusion Energy 2008[END_REF], NIMROD [START_REF] Izzo | [END_REF] [19] [20] [21] [22], JOREK [23] and the code of Rozhansky et al. [24]. Compared to these works, the approach presented here is original in that it treats gas transport as purely convective (as it should be, according to first principles) and it includes the gas reservoir, vacuum region and plasma in the simulation domain. These ingredients are essential to simulate the processes described below.

The paper is constructed as follows. Section II presents the model and the IMAGINE code. Section III describes simulations of a D 2 MGI into a JET Ohmic plasma, discusses the mechanisms at play and compares the results to experimental observations. Section IV presents simulations whose aim is to investigate the penetration of the gas into the background plasma co-existing with a RE beam. Finally, Section V concludes and gives perspectives.

II. MODEL

A. Geometry, equations and assumptions

The model is 1D in the radial direction and assumes a slab geometry. It is a fluid model, which seems reasonable given the large gas densities typical of MGI. The gas is assumed to be made of atoms of the same species as the plasma ions, with an atomic number Z = 1. Below, we simulate D 2 MGI into D + plasmas. We however ignore the existence of D 2 molecules and do as if the gas was made of D atoms, neglecting the dissociation energy of D 2 molecules (this seems reasonable since this energy is a factor 5 smaller than the ionization energy of two D atoms). The model consists of 6 equations describing the evolution of the electron density n e (which, due to quasi-neutrality, is equal to the ion density n i ) and temperature T e , the ion temperature T i and the neutral density n n , velocity V n and pressure P n , written here in conservative form:

∂ t n e = n e n n I -n 2 e R + ∂ r (D∂ r n e ) (1)
∂ t ( 3 2 n e eT e ) = -n e (n n IE ion + n n L lines + n e R 3 2 eT e ) -n 2 e L brem+rec + ∂ r (χn e ∂ r (eT e )) ( 2
)
∂ t ( 3 2 n e eT i ) = 3 2 n e (IP n -n e ReT i -σ cx V cx (n n eT i -P n )) + ∂ r (χn e ∂ r (eT i )) (3)
∂ t n n = -∂ r (n n V n ) -n e n n I + n 2 e R (4)
∂ t (m n n n V n ) = -∂ r (m n n n V 2 n + P n ) -n n n e (I + σ cx V cx )m n V n (5)
∂ t ( 3 2 P n + 1 2 m n n n V 2 n ) = -∂ r ( 5 2 P n V n + 1 2 m n n n V 3 n ) -n n n e (I + σ cx V cx )( 3 2 P n /n n + 1 2 m n V 2 n ) + n e (n e R + n n σ cx V cx ) 3 2 eT i (6)
In these equations, I, R and σ cx V cx are the ionization, recombination and charge exchange rates, E ion is the ionization energy and L lines and L brem+rec are the line and continuum (i.e.

bremmstrahlung plus recombination) radiation rates. Details on these coefficients are given in Section II C. All quantities are in SI units except temperatures which are in eV. D and χ are particle and heat diffusivities which are meant to represent turbulent transport. These two parameters are the only ad hoc parameters of the model and their effect is in fact limited.

In the simulations below, we use

D = χ = 1 m 2 /s.
It is worth comparing our equations to those of the general non-linear fluid model for reacting plasma-neutral mixtures from Meier and Shumlak [26]. Equations 1-6 above correspond to Equations 48-56 of [26] under the following assumptions/modifications: charged species are assumed to be at rest; scattering collisions are neglected; radiation losses are added; the term R cx in is neglected; electron particle diffusion is added and electron and ion heat fluxes are treated as diffusive as a way to account for turbulent transport, as mentioned above.

The assumption of charged species being at rest is probably the most important one and deserves some discussion. As mentioned by Rozhansky et al. [24], the gas tends to "push" charged species, i.e. to give rise to an E × B flow in the direction of the gas flow, but several mechanisms may damp this E × B flow. The physics involved is similar to the physics of pellet ablation clouds [25]. In the latter case, the strongest damping mechanism may be described as follows. The E × B flow is intimately related to an electric field E = -∇Φ in the cloud. If a field line connects the top and the bottom of the cloud (after a certain number of toroidal turns), the potential difference ∆Φ between its extremities will drive a parallel current. The current path will be closed (it has to be for quasi-neutrality to pertain) by a polarization current proportional to dE/dt circulating inside the cloud. The latter is associated to a reduction of |E| and therefore a damping of the flow. This mechanism, which is already strong for pellet ablation clouds, is likely to be even stronger in the MGI case because of the larger the cloud poloidal extent. Note that this mechanism, introduced by Pégourié et al. in [25], was not included in the slightly older work of Rozhansky et al. [24].

For simplicity we assume here that this braking mechanism is strong enough that charged species remain at rest. The question may however be addressed in future work by including a momentum equation for charged species in our model. Looking at our equations, one may notice that the term -n e (I+ < σv > cx ) 1 2 m n n n V 2 n in the neutral energy equation (Eq. 6), which represents kinetic energy transferred from neutrals to ions due to ionization and charge exchange, has no counterpart in the ion energy equation (Eq. 3). This is because this term is associated to a momentum transfer between neutrals and ions, and this momentum is assumed to be annihilated by the braking force(s) acting on charged species. Consistently, the associated energy is assumed to be annihilated by the work of the braking force(s).

B. Initial and boundary conditions

One key feature of our model is to treat in a unified way the plasma, the vacuum region and the gas reservoir: the above equations apply throughout these three domains and the simulations cover the gas propagation into the vacuum and its interaction with the plasma, as will be presented in Section III.

Boundary conditions are null gradients for all fields at the center of the plasma (r = 0) and V n = 0 at the end of the gas reservoir (r = r max ).

Initial conditions for n e , T e and T i are typically set according to pre-MGI experimental measurements. Initial conditions for n n and P n are a finite and homogeneous value in the reservoir and 0 elsewhere, while V n = 0 everywhere. The initial values of n n and P n in the reservoir, as well as the radial extent of the reservoir in the simulation, are chosen so as to match three critical experimental quantities: the initial number of D atoms and the sound velocity c s,res inside the reservoir and the initial flux of atoms out of it (technical details are given in Appendix A). This allows reproducing the gas flow in the vacuum region, as we shall see below.

C. Atomic physics coefficients

Coefficients I, R, L lines and L brem+rec are taken from the ADAS database [START_REF]Atomic Data and Analysis Structure User Manual[END_REF], where they are given as function of n e and T e . They are represented in Figure 1.

FIG. 1. Atomic physics coefficients for deuterium from the ADAS database

For the charge exchange rate σ cx V cx , we use the expressions given in Meier and Shumlack [26] (which apply to hydrogen but we assume that deuterium is similar):

σ cx [m 2 ] = 1.09 × 10 -18 -7.15 × 10 -20 ln(V cx ) (see Section III.D.2 in [26]) and V cx = ( 4 π v 2 T i + 4 π v 2 T n + V 2 n ) 1/2 (
Eq. 17 in [26]), where v T α (2eT α /m α ) 1/2 is the thermal velocity of species α. Figure 2 shows σ cx V cx as a function of T i , assuming that v T n and V n are negligible. the neutral equations ignoring the atomic physics terms. The MUSCL scheme is a finite volume method providing highly accurate numerical solutions even when the solutions exhibit shocks or discontinuities (which is the case here). The part related to atomic physics is advanced separately with an explicit scheme. The equations for electrons and ions are advanced with an explicit conservative scheme.

III. SIMULATION OF A D 2 MGI INTO A JET OHMIC PLASMA

In this Section, we present IMAGINE simulations of MGI into a "hot" plasma. Section IV will deal with simulations of MGI into a "cold" plasma.

A. Description of the experiment

We simulate JET pulse 86887. This is an Ohmic D plasma pulse with B t = 2 T, I p = 2 MA, q 95 = 2.9 in which a disruption was triggered on a "healthy" plasma by activating the Disruption Mitigation Valve number 2 (DMV2). Note that the same pulse has been modelled with the 3D non-linear MHD code JOREK, as reported in [23].

B. Simulation results

We shall now describe the IMAGINE simulation results, starting with the propagation of the gas into the vacuum region (i.e. the guiding tube).

Gas propagation into the vacuum region

In laboratory experiments, Bozhenkov et al. have found that the gas flow in the guiding tube is well described by the analytical solution of the 1D Euler equations with an infinite reservoir, under the condition that d cs,res•t > 0.6, where d is the distance to the reservoir [9]. This type of flow is called a "rarefaction wave". Since in the vacuum region the equations of IMAGINE boil down to the 1D Euler equations, we recover this flow in the simulations.

Figure 6 (left) shows successive neutral density profiles calculated by IMAGINE. One can observe the typical self-similar behaviour of the rarefaction wave [9]. Figure 6 (right) shows the flux of D 2 molecules across r = 3 m, comparing IMAGINE and the analytical solution (given explicitly in [9]). A good match (with a difference < 10 %) can be observed up to about 3 ms, from which point the IMAGINE solution starts to decrease due to the depletion of the reservoir, in contrast with the analytical solution which assumes an infinite reservoir.

In Figure 6 (left), the reader may notice that the reservoir in IMAGINE has a radial extent of about 1.5 m, which is much larger than the typical dimension of the actual reservoir, whose volume is V res = 1 L. Reasons for this are detailed in Appendix A. One consequence is that the pressure equilibration time in the reservoir is much longer in the simulations than in the experiment. This must lead to a slightly inexact description of the gas flow. The real gas flux probably has a smoother behaviour, with no abrupt change at 3 ms as in Figure 6 (right). However, the discrepancy is probably not very large and, as stated above, the early gas flow calculated by IMAGINE matches well the laboratory measurements of Bozhenkov et al., giving confidence in the following results (which essentially have to do with the early interaction of the gas with the plasma). Let us now focus on what happens when the gas reaches the plasma. Here, charge exchange and recombination play a central role. In order to demonstrate this, we will first present a simulation without these effects and then a simulation with these effects.

2. Gas penetration into the plasma, neglecting charge exchange and recombination

Without charge exchange and recombination, the gas penetrates easily into the plasma. This is apparent in Figure 7 (top left), which shows profiles of the neutral velocity V n at different times: the gas flow is not affected as the gas enters the plasma. On the neutral density profiles shown in Figure 7 (top right), one can see that the tip of the rarefaction wave is "eaten" by ionisation. This causes a dramatic increase in n e , which reaches values on the order of 2 • 10 21 m -3 , as can be seen in Figure 7 (bottom left) (note the logarithmic scale in this figure). This increase in n e is clearly incompatible with experimental data, as can be seen in Figure 8, which compares experimental and synthetic interferometry data. In this figure, the synthetic data is calculated from the output of IMAGINE as n el = 2 a 0 n e dr (the factor 2 comes from the fact that the experimental measurement corresponds to a -a n e dr = 2 a 0 n e dr), and compared with experimental measurements from lines of sight 2 and 3 of the interferometer, which are rather central, as can be seen in Figure 8 (right). Note that the comparison is meaningful only if the experimental n e is relatively homogeneous on flux surfaces. This may indeed be the case due to the fast parallel expansion of the over-density created by the MGI. For example, JOREK simulations [23] find a parallel expansion velocity v 30km.s -1 , corresponding to a density homogenization time τ 2πqR v 2ms. In Figure 8, the synthetic signal from the IMAGINE simulation without charge exchange and recombination goes beyond the axis limit at about 2.5 ms and continues increasing dramatically later, reaching values one order of magnitude higher than the experimental ones, clearly indicating that some assumptions of the model must be wrong.

This process is accompanied by the penetration of a very sharp cold front, as can be seen in Figure 7 (bottom right). The cooling is mostly due to dilution, since n e increases by roughly 2 orders of magnitude. It can be seen that within 1.5 ms from the DMV2 trigger, the cold front has reached r = 0.9 m, which corresponds to the q = 2 surface, and within 4 ms it has reached mid-radius. This is also incompatible with experimental data since the TQ is observed about 12 ms after the DMV2 trigger. Recall that in several machines, the TQ has been observed to occur when the cold front reaches the q = 2 surface [6] [8] [7], and even though there is not clear experimental evidence for this fact in JET, JOREK modelling suggests that this is the case too [23]. To summarize, when neglecting charge exchange and recombination, the gas and the cold front penetrate much too fast and easily compared to the experiment.

Gas penetration into the plasma, including charge exchange and recombination

Turning on the charge exchange and recombination terms changes the situation drastically: the gas flow is now rapidly stopped at the edge of the plasma and a shock wave is created which propagates away from the plasma, as is clearly visible on the successive neutral velocity (V n ) and density (n n ) profiles from Figure 9 (top left and right plots). Concerning the shock wave, it is good to check whether Rankine-Hugoniot jump conditions [START_REF] Landau | [END_REF] are satisfied in the simulation. Figure 10 shows that this is the case: while n n is discontinuous across the shock (top left plot), the quantities This indicates that IMAGINE resolves the shock wave appropriately.

ρ n (V n -V shock ) (top right), ρ n V n (V n -V shock ) + P n (bottom left) and (V n -V shock )( Pn γ-1 + 1 2 ρ n V 2 n ) + V n P n (
As a result of the strong braking, only a small fraction of the gas penetrates into the plasma and the n e increase is therefore much more modest than in the previous simulation.

Still, n e increases by a factor 6 at r = 0.9 m, as can be seen in Figure 9 (bottom left plot). Looking at Figure 8 again, this simulation (plain line) appears much more consistent with interferometry data than the previous one: the n el increase now has the right order of magnitude. Figure 9 (bottom right) indicates that the cold front is less sharp and penetrates much more slowly than in the previous simulation. It takes about 5 -10 ms for it to reach the q = 2 surface (r = 0.9 m), which is compatible with the experimental TQ onset time of 11 ms. In conclusion, when including charge exchange and recombination, IMAGINE simulations are rather consistent with experimental measurements. Let us now interprete these results.

C. Interpretation

The above simulations demonstrate the crucial role of charge exchange and recombination in promptly braking the gas, limiting its penetration into the plasma and giving birth to a shock wave propagating away from the plasma. In order to identify the mechanisms responsible for this braking, let us consider the equation for the evolution of V n (which can be derived from Equations 4 and 5):

m n n n d t V n = -∂ r P n -(n 2 e R + n e n n σ cx V cx )m n V n (7)
Here,

d t = ∂ t + V n ∂ r .
This equation shows that the neutral velocity is influenced by

(V n -V shock ) (top right), ρ n V n (V n -V shock ) + P n (bottom left) or (V n -V shock )( Pn γ-1 + 1 2 ρ n V 2 n ) + V n P n (bottom right),
in agreement with the RH conditions.

the neutral pressure gradient and by a friction force (n 2 e R + n e n n σ cx V cx )m n V n . The latter comes from the fact that each charge exchange or recombination creates a neutral which, on average, is at rest, thereby decreasing the average neutral velocity V n by a quantity proportional to V n /n n . Regarding the former, let us consider the equation for the evolution of P n (which can be derived from Equations 4, 5 and 6):

∂ t P n = -V n ∂ r P n - 5 3 P n ∂ r V n + n n n e σ cx V cx ( 1 3 m n V 2 n + eT i - P n n n ) + n 2 e R( 1 3 m n V 2 n + eT i ) -n e IP n (8)
It appears that charge exchange and recombination tend to increase P n via (among others) terms proportional to eT i . This is due to the fact that each charge exchange or recombination creates a neutral whose energy, on average, is equal to 3 2 eT i . This effect typically leads to the formation of a P n spike at the gas-plasma interface. The -∂ r P n term in Equation 7then tends to accelerate the gas on the inside of this spike and decelerate it on the outside.

From a kinetic point of view, this may be understood as hot neutrals being created with a velocity pointing randomly either towards the plasma core or away from it (depending on the gyro-angle of the ion just before the charge exchange or recombination).

In summary, charge exchange and recombination influence the gas flow in two ways: via a friction force and via gas heating. The relative importance of these two mechanisms depends on the ion temperature. We will come back to this point in Section IV.

The fact that a shock wave should appear when the gas reaches the plasma has been predicted already by Parks and Wu [29], who state (at the beginning of Section 3 of [29]) that "The plasma pressure impresses on the frontal surface driving a shock wave running inwards and backwards". In our opinion, this statement may be misleading since the plasma pressure does not directly apply on the gas. However, the gas heating effect described above may be seen as a conversion of plasma pressure into neutral pressure and in that sense, plasma pressure may be considered to indirectly impress on the gas. But the above discussion shows that this is not the only type of interaction between the plasma and the gas: the frictional interaction should not be forgotten.

Another question is whether there is a dominant player between charge exchange and recombination. A partial answer is provided in Figure 8 (left), which contains simulation results with charge exchange and recombination turned off separately. This figure indicates that charge exchange has a stronger effect than recombination, but the most important point is probably that each of these two effects is strong enough to significantly brake the gas and create a shock wave by itself.

IV. GAS PENETRATION INTO THE BACKGROUND PLASMA SURROUND-ING A RUNAWAY BEAM

As mentioned in the introduction, multiple attempts to "kill" an RE beam with MGI (including with high Z gases like Kr or Xe) in JET have shown no effect [1]. A possible explanation may be that the RE beam, which is localized in the core of the plasma, extending perhaps up to mid-radius, is "shielded" by the cold background plasma surrounding it via a gas braking effect similar to the one described above. However, it is essential to note that on Tore Supra [30], DIII-D [31] or ASDEX-U [START_REF] Pautasso | Generation and suppression of runaway electrons in ASDEX Upgrade disruptions[END_REF], MGI clearly affects the RE beam, so the gas must reach the beam in these machines. It is therefore interesting to investigate possible reasons for differences in terms of gas penetration between the machines.

The electron density n e,bg and temperature T e,bg of the background plasma are one possibility. Although measurements are not easy (especially for temperature), n e,bg seems to be typically a few 10 20 m -3 in JET and DIII-D and a few 10 19 m -3 in Tore Supra and ASDEX-U, and T e,bg is estimated to be at least 20 eV in JET while T e,bg ∼ 2eV in DIII-D and ASDEX-U [1] [31] [START_REF] Pautasso | Generation and suppression of runaway electrons in ASDEX Upgrade disruptions[END_REF]. The background plasma therefore seems to be denser and/or hotter in JET than in DIII-D or ASDEX-U. IMAGINE allows investigating how gas penetration depends on n e,bg and T e,bg . For this purpose, we made simulations of D 2 MGI in JET from DMV2 at 20 bar. This setting was chosen because it corresponds to an actual pulse. The n e and T e profiles were scaled from the profiles of the Ohmic pulse 86887 modelled above. We assumed T i = T e for the initial profiles. We also turned off the radiation and recombination terms. Indeed, due to the low T e and to the absence of certain physical effects in the model (e.g. Joule heating), the bulk plasma may otherwise undergo a radiative collapse and/or significant recombination over the timescale of the simulation, which are not the purpose of the present study. Figure 11 Simulation (A) shows that even at very low T e,bg , the gas is prevented from penetrating far into the plasma if n e,bg is large enough. The responsible mechanism for this is gas-plasma friction due to charge exchange, which remains significant even at very low plasma temperatures due to the weak T i dependency of σ cx V cx at low T i , visible in Figure 2 (right). At lower densities, however, such as in simulation (B), the gas may penetrate much further. This is easily understood as the friction force is proportional to n e . Simulations (B), (C) and (D) constitute a scan in T 0 e,bg at fixed n 0 e,bg = 10 19 m -3 . Moving from T 0 e,bg = 2 eV to T 0 e,bg = 20 eV, i.e. from (B) and (C), has almost no influence on the result, whereas at T 0 e,bg = 2 keV (case (D)), gas penetration is clearly hindered. This is due to the T i dependency of the above-mentioned mechanism of gas braking via charge exchange heating: for T 0 e,bg < 20 eV, this mechanism plays an almost negligible role, whereas for T 0 e,bg = 2 keV, it is the dominant player. The simulations therefore reveal an important difference in gas penetration physics between a hot and a cold target plasma: in the former case, the gas is mainly braked via heating while in the latter, it is mainly braked by friction.

It may be that the situation at JET corresponds to simulation (A), where the superficial penetration of the gas would be insufficient to affect the RE beam if the latter is localized inside, say, r/a 0.5. The lower n 0 e,bg at Tore Supra and ASDEX-U may on the other hand correspond to simulation (C), where one would expect the gas to affect the beam. Simulation results however do not provide an obvious explanation for the difference between JET and DIII-D, since n 0 e,bg has the same order of magnitude in these two machines and the factor ∼ 10 difference in terms of T 0 e,bg is not expected to affect gas penetration significantly. Of course, it may be that the explanation lies somewhere else, for example in a machine size effect (a larger machine means a thicker background plasma "shield" around the RE beam) or different MGI settings. Investigating these possibilities is left for future work. For now, the important practical conclusion of our study in view of further experiments in JET is that optimizing the RE beam production "recipe" in order to lower n 0 e,bg may help get an effect of MGI on the beam.

To finish this section, we would like to draw attention to an interesting observation made in the above-mentioned Tore Supra experiment. Figure 7b in [30] shows that the neutral pressure measured with a pressure gauge in the vacuum vessel decreases abruptly is confined to a smaller volume and/or hotter than 120 • C. This is consistent with the findings presented in this paper that gas-plasma interaction prevents most of the gas from penetrating the plasma volume (even when part of the gas does penetrate the plasma and affects the runaway beam) and also heats up the gas. An interesting direction for future work is to investigate the possible role of this gas surrounding the plasma. This question cannot be addressed by a simple 1D model like the one used in IMAGINE.

V. CONCLUSION AND PERSPECTIVES

The model introduced in this paper is simple in several respects: it uses a 1D radial slab geometry, treats only D 2 MGI into a D + plasma, and assumes ionized species to be at rest.

However, the model comprises features which are essential to investigate the fundamental physics of gas penetration during an MGI: the gas equations are based on first principles, and the simulation domain encompasses the plasma, vacuum region and gas reservoir.

IMAGINE simulations have revealed that the gas delivered by MGI is in general strongly braked as it reaches the plasma. This occurs due to two different mechanisms: a gas-plasma friction force and a heating of the gas by the plasma. Both mechanisms are related to charge exchange and recombination. The heating mechanism typically dominates in high temperature, pre-TQ plasmas, whereas the friction mechanism dominates in cold, post-TQ plasmas. The strong gas braking results in a partial penetration of the gas into the plasma and provokes the appearance of a shock wave in the gas, which propagates away from the plasma, braking and compressing the incoming gas. It is important to note that, although the IMAGINE model is presently restricted to D 2 MGI into a D + plasma, we expect these mechanisms to generalize to cases where the gas is of a different species than the plasma ions. Implementing other gases than D 2 in MGI is a direction for future work. Another important question for future research is whether the hypothesis of ionized species being at rest is justified or whether the gas may "push" the plasma.

Simulation results are quantitatively consistent, at least in terms of orders of magnitude, with experimental observations for a D 2 MGI in an Ohmic JET pulse. We stress that no parameter was adjusted in order to match experimental measurements. In fact, the model does not contain adjustable parameters besides the diffusion coefficients D and χ, which only have a small effect on the results. More simulations and comparisons to experiments are necessary to assess whether the model, in spite of its simple geometry, really has a predictive capability or whether the agreement reported in this paper is fortuitous.

IMAGINE allows investigating possible differences between machines regarding experiments aiming at RE beam mitigation with MGI. Simulations show that the gas penetration depth is a decreasing function of the electron density of the target plasma. It may be that at JET, the density of the background plasma co-existing with the RE beam, n 0 e,bg , is too high for the gas to be able to reach the beam, which would explain its absence of effect. In contrast, in Tore Supra and ASDEX-U, n 0 e,bg is one order of magnitude lower and the gas may penetrate, consistently with the positive experimental results. However, in DIII-D, in spite of a value of n 0 e,bg comparable to JET, a positive result is obtained. It may be noticed that T 0 e,bg is one order of magnitude lower in DIII-D than in JET, but IMAGINE simula-tions indicate that this should not affect gas penetration much. Further investigations are therefore needed to identify the origin of the difference between JET and DIII-D.

Another important area for future work is to implement the physics treated in IMAGINE in other models used for MGI simulations, e.g. 3D non-linear MHD codes such as JOREK or NIMROD: the gas deposition model is indeed rather ad hoc at present, which limits their predictive capability. It may be surprising that the reservoir pressure in the simulation is so low compared to the experimental one (2.48 Pa vs. 5 bar). However, it should be kept in mind that in our model, the plasma is as completely surrounded by the gas reservoir and the incoming gas flow is spread over the whole plasma surface instead of being very localized.

FIG. 2 .

 2 FIG. 2. D-D+ charge exchange rate σ cx V cx as a function of T i (left: linear scale; right: logarithmic scale), assuming that v T n and V n are negligible (see text for explanation)

FIG. 3 .FIG. 5 .

 35 FIG. 3.Experimental n e (left) and T e (right) profiles for JET pulse 86887 from High Resolution Thomson Scattering (HRTS) just before the DMV2 trigger (red crosses) and fits of these profiles used as initial conditions in IMAGINE (blue). (Note: the HRTS data has been remapped onto an EFIT equilibrium. The radial coordinate r used in the above profiles is defined as a • ψ 0.5 n , with a = 1 m, where ψ n is the normalized poloidal flux from EFIT.)

FIG. 6 .

 6 FIG. 6. Left: Successive neutral density profiles (separated by 0.1 ms) from IMAGINE simulation. Right: Flux of D 2 molecules in vacuum across r = 3 m as found with IMAGINE and with the analytical solution of the 1D Euler equations.

FIG. 7 .

 7 FIG. 7. Successive profiles of neutral velocity V n (top left), neutral density n n (top right), and electron density n e (bottom left), separated by 0.5 ms; electron temperature T e versus time and radius (bottom right)

FIG. 8 .

 8 FIG.8. Comparison between experimental and synthetic interferometry data (see Figure4to visualize the interferometry lines of sight).

  bottom right) are continuous (where ρ n = m n n n ), in agreement with the conservation of mass, momentum and energy density.

FIG. 9 .

 9 FIG. 9. Successive profiles of neutral velocity V n (top left), neutral density n n (top right), and electron density n e (bottom left), separated by 0.5 ms; electron temperature T e versus time and radius (bottom right)

FIG. 10 .

 10 FIG. 10. Verification of the Rankine-Hugoniot (RH) jump conditions: the shock (blue dashed line) corresponds to a discontinuity on the neutral density (top left) but not on ρ n (V n -V shock ) (top

  presents the evolution of the sum of free plus bound electrons, i.e. n e + n n , versus time and radius, for four simulations: (A) n 0 e,bg = 10 20 m -3 , T 0 e,bg = 2 eV (top left); (B) n 0 e,bg = 10 19 m -3 , T 0 e,bg = 2 eV (top right); (C) n 0 e,bg = 10 19 m -3 , T 0 e,bg = 20 eV (bottom left); (D) n 0 e,bg = 10 19 m -3 , T 0 e,bg = 2 keV (bottom right). The superscript 0 denotes the central value at the beginning of the simulation. We chose to represent n e + n n because both free and bound electrons may brake RE (although not exactly to the same extent).

from 2 . 5 FIG. 11 .

 2511 FIG. 11. Electron plus neutral density as a function of time and radius for n 0 e,bg = 10 20 m -3 , T 0 e,bg = 2 eV (top left); n 0 e,bg = 10 19 m -3 , T 0 e,bg = 2 eV (top right); n 0 e,bg = 10 19 m -3 , T 0 e,bg = 20 eV (bottom left); n 0 e,bg = 10 19 m -3 , T 0 e,bg = 2 keV (bottom right)

n sim,D res = 1 .

 1 43 • 10 21 m -3 , P sim,D res = 2.48 Pa and δ sim r,res = 1.43 m.

VI. ACKNOWLEDGEMENTS

We would like to thank H. Bufferand, E. Joffrin, B. Pégourié and C. Sommariva for helpful discussions or comments on the manuscript. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the ITER Organization or the European Commission.

Appendix A: Details on simulation settings

In this Appendix, we detail how the initial density n sim res and pressure P sim res and of the reservoir and its radial extent δ sim r,res are set so as to match the initial number of D atoms in the reservoir N res , the initial flux of D atoms out of it Φ res , and the initial sound velocity in the reservoir c s,res .

One should keep in mond that in the simulations, n n is the D atom density but in reality the gas is made of D 2 molecules. In the equations below, we use a superscript to specify whether quantities refer to D atoms or D 2 molecules.

Considering the simulation domain as a slab of length 2πR 0 and height 2πa, with R 0 the major radius and a the minor radius of the machine, the three above-mentioned conditions translate to: Equation A2 anticipates on the result presented in the next section that the gas velocity at the exit of the DMV is the sound velocity (both in the experiment and simulation).

After some simple algebra, Equations A1, A2 and A3 yield the following expressions for the input parameters n sim,D res , P sim,D res and δ sim r,res :

where V exp res is the volume of the gas reservoir in the experiment. In JET pulse 86887, the experimental parameters are P exp,D


E Nardon

A Fil

M Hoelzl

G Huijsmans

Progress in understanding disruptions triggered by massive gas injection via 3D non-linear MHD modelling with JOREK

3D non-linear MHD simulations of a D 2 Massive Gas Injection (MGI) triggered disruption in JET with the JOREK code provide results which are qualitatively consistent with experimental observations and shed light on the physics at play. In particular, it is observed that the gas destabilizes a large m/n = 2/1 tearing mode, with the island O-point coinciding with the gas deposition region, by enhancing the plasma resistivity via cooling. When the 2/1 island gets so large that its inner side reaches the q = 3/2 surface, a 3/2 tearing mode grows. Simulations suggest that this is due to a steepening of the current profile right inside q = 3/2. Magnetic field stochastization over a large fraction of the minor radius as well as the growth of higher n modes ensue rapidly, leading to the Thermal Quench (TQ). The role of the 1/1 internal kink mode is discussed. An Ip spike at the TQ is obtained in the simulations but with a smaller amplitude than in the experiment. Possible reasons are discussed.

Introduction

In order to develop a predictive modelling capability in support of the design and operation of the ITER Disruption Mitigation System (DMS) [1] [2], it is essential to validate models on present devices as well as to understand the physical mechanisms at play. The present paper describes progress made with the JOREK 3D non-linear magnetohydrodynamics (MHD) code in the domain of Massive Gas Injection (MGI), one of the options considered for the ITER DMS. The focus is put on the pre-Thermal Quench (TQ) and TQ phases. The Current Quench (CQ) phase will be the object of future work.

A publication by Fil et al. in 2015 [3] introduced first results from simulations of a pure D 2 MGI into a JET Ohmic plasma. Work has been pursued since then on the same case and with essentially the same model, but with more realistic input parameters thanks to numerical improvements. This has lead to significant progress which is presented in this paper. Most notably, while only an incomplete TQ was obtained in [3], simulations now display a more proper TQ. Furthermore, numerical experiments as well as a detailed analysis of simulation results have led to a better understanding of the physics at play.

The paper is constructed as follows. In Section 2, the experiment is introduced. Section 3 presents the model and Section 4 the setup of simulation parameters, highlighting differences with [3]. Section 5 gives an overview of a simulation, after which Sections 6 and 7 describe investigations on the mechanisms at play during the pre-TQ and TQ phases, respectively. Finally, Section 8 concludes.

The experiment

The modelled case has already been described in the paper by Fil et al. [3] but it is briefly described again here for convenience. This is JET pulse 86887, an Ohmic pulse with B t = 2 T, I p = 2 MA, q 95 = 2.9 in which a disruption was MGI-triggered on a "healthy" plasma by activating the Disruption Mitigation Valve number 2 (DMV2) pre-loaded with D 2 at 5 bar. The gas is injected at the outer midplane. Central values of the electron density and temperature prior to the MGI are n e = 3 • 10 19 m -3 and T e = 1.2 keV. The volume of the DMV2 reservoir is 10 -3 m 3 and its temperature is about 300 K, so it initially contains about 1.2 • 10 23 D 2 molecules, which represents roughly 100 times more D nuclei than initially present in the plasma.

Figure 1 shows an overview of the disruption phase. First effects of the MGI are visible from about 2 ms (relative to the DMV2 trigger) in the form of increases in the line integrated density and radiated power. The TQ occurs at about 12 ms as can be seen from the fast collapse of the central Soft X-Ray (SXR) signal accompanied by a burst of MHD activity and immediately followed by the characteristic I p spike. The CQ ensues. Experimental time traces for JET pulse 86887, from top to bottom: plasma current Ip, magnetic fluctuations from Mirnov coil, radiated power from bolometry, line integrated density from interferometry, and soft X-rays signal from a central chord. The time origin corresponds to the DMV2 trigger.

The model

The model used here is essentially the same as in [3] (with only a few differences detailed below). The reader is therefore referred to [3] for a more complete model description. Nevertheless, for convenience we write again the equations here:

∂ψ ∂t = η∆ * ψ -R [u, ψ] -F 0 ∂u ∂φ (1)
j = ∆ * ψ (2) R∇ • R 2 ρ∇ pol ∂u ∂t = 1 2 R 2 |∇ pol u| 2 , R 2 ρ + R 4 ρω, u + [ψ, j] - F 0 R ∂j ∂φ + ρT, R 2 + Rµ∇ 2 ω +∇ • ρρ n S ion -ρ 2 α rec R 2 ∇ pol u (3) ω = ∇ 2 pol u = 1 R d dR (R du dR ) + d 2 u dZ 2 (4)
∂ρ ∂t = -∇•(ρv)+∇•(D ⊥ ∇ ⊥ ρ+D ∇ ρ)+ρρ n S ion -ρ 2 α rec (5)
Progress in understanding disruptions triggered by massive gas injection via 3D non-linear MHD modelling with JOREK3

∂(ρT ) ∂t = -v • ∇(ρT ) -γρT ∇ • v +∇ • (κ ⊥ ∇ ⊥ T + κ ∇ T ) + 2 3R 2 ηj 2 -ξ ion ρρ n S ion (6)
ρB 2 ∂v ∂t = -ρ F 0 2R 2 ∂(B 2 v 2 ) ∂φ - ρ 2R B 2 v 2 , ψ - F 0 R 2 ∂(ρT ) ∂φ + 1 R [ψ, ρT ] + B 2 µ (T )∇ 2 pol v +(ρ 2 α rec -ρρ n S ion )B 2 v ( 7
)
∂ρ n ∂t = ∇ • (D n : ∇ρ n ) -ρρ n S ion + ρ 2 α rec + S n (8)
ψ is the poloidal flux, j = Rj φ the toroidal current density times the major radius, u the poloidal flow potential, ω the toroidal vorticity, ρ the plasma mass density, T the total (ion + electron) temperature, v the parallel (to the magnetic field) velocity and ρ n the neutral mass density. S ion and α rec designate respectively the ionization and recombination rate coefficients for deuterium, parameterized as detailed in [3]. Note that apart from the neutral source, there are no sources in the model. Realistic sources indeed have a negligible effect on the small timescale of the disruption.

Hyper-resistivity/conductivity/viscosity terms, not shown above, are also included for numerical stability reasons.

There are two differences between the present model and that of [3]. First, for the sake of energy conservation, the resistivity used in the Ohmic heating term of Eq. 6 is the same as in Eq. 1, while in [3] the Spitzer resistivity was used. Second, D radiation terms have been suppressed for simplicity after verifying that they do not influence the results significantly.

Simulations setup

Initial conditions and general parameters

Initial density and temperature profiles are taken from Thomson scattering measurements, asssuming T i = T e . The initial magnetic equilibrium is calculated from EFIT data as explained in [3]. Toroidal rotation is initialized at 0. Table 1 gives the value of input parameters used here as well as in [3]. Improvements in the poloidal mesh now allow running simulations with a hyper-resistivity η h and hyper-viscosity µ h roughly 500 times smaller than in [3], which removes the artificial stabilization of the MHD activity due to these coefficients in [3]. Also, present simulations use much smaller diffusion coefficients than in [3], both for ions/electrons and for neutrals. The artificial core dillution observed in [3] is therefore not present anymore. It is important to recall that realistic temperature dependencies are used for the resistivity and parallel thermal conductivity: η = η 0 • (T 0 /T ) 3/2 and κ = κ 0 • (T /T 0 ) 5/2 , where T 0 is the initial temperature at the center of the plasma. For κ 0 , the Spitzer-Härm value is used, while η 0 is typically one order of magnitude larger than the Spitzer value.

Gas source setting

The physics of gas penetration into the plasma during an MGI is rather complex [4] and it is clear that the present model is not appropriate to describe it accurately. In other terms, JOREK simulations are not predictive in this respect. The gas source is therefore adjusted so as to best match interferometry data. Its temporal shape is based on the solution for the gas flow in vacuum [3] but using a DMV2 pressure of 0.1 bar instead of the experimental 5 bar (which can be interpreted as a sign of the low fuelling efficiency of the MGI). Note that in [3], a much larger value of typically 1 bar was needed due to the use of much larger particle diffusion coefficients. The gas source is located at the edge of the plasma (centered 5 cm inside the separatrix) and localised both poloidally and toroidally (see [3] for precisions). Figure 2 shows experimental and synthetic line integrated densities n el for three interferometry Lines of Sight (LoS). Three simulations are shown, with different toroidal localisations of the source: ∆φ M GI = 0.5 (dashed), 2 (plain) and 100 (i.e. virtually axisymmetric, dashdotted) respectively (angles are given in radians). Even though the same number of neutrals is injected in each simulation, clear differences appear in the synthetic n el . In order to interprete these differences, one should consider that interferometry chords are located 180 degrees away toroidally from DMV2, as shown in Fig. 3. Poloidal cross-sections of n e in the interferometer plane are shown in Fig. 4 for ∆φ M GI = 0.5 (top) and ∆φ M GI = 100 (bottom). These figures indicate that LoS 2 (red) runs through a region which is connected to the gas deposition region and therefore sees a faster increase in n el than LoS 3 (blue) for ∆φ M GI = 0.5 (dashed) but not for ∆φ M GI = 100 (dash-dotted). The simulation at ∆φ M GI = 0.5 actually matches experimental data very well for the three LoS until about 3 ms, after which the agreement deteriorates for a yet unidentified reason. Over a larger period of time, the case ∆φ M GI = 2 is the one that offers the best global match and is therefore selected for a detailed description and analysis below. mode (q 0 < 1 in this simulation, consistently with the Progress in understanding disruptions triggered by massive gas injection via 3D non-linear MHD modelling with JOREK5

presence of sawteeth in this pulse). In the third row, at t = 5.7 ms, the magnetic field has become stochastic over roughly the outer half of the plasma, which leads to a temperature flattening by parallel conduction in this region. However, the core is still hot at this stage. Finally, at t = 6.2 ms (fourth row), the core T e has collapsed via a convective mixing related to the 1/1 mode. The behaviour observed in this simulation is similar to the one found in NIMROD simulations for other tokamaks [6] [7].

Physics of the pre-thermal quench phase

Let us now analyse the mechanisms at play, beginning with the pre-TQ phase. An important question is how the gas generates the 2/1 tearing mode. One may imagine at least three possible mechanisms. The first one is that the gas a 3D pressure field to which j and B have to adjust in order for force balance to pertain, leading to a 3D equilibrium which may imply a 2/1 tearing mode. The other two mechanisms, in contrast to the first one, are directly related to the temperature dependence of the resistivity η. One is an axisymmetric effect: the gas cools the edge of the plasma, increasing η and therefore contracting the current profile. This should increase the current density gradient |j φ | at the edge of the still hot region. When the steep gradient gets just inside the q = 2 surface, the 2/1 tearing mode should be destabilized via the following term in the resistive MHD energy principle [START_REF] Biskamp | Nonlinear Magnetohydrodynamics[END_REF]: δW ∝ j φ 1-nq/m |ψ 1 | 2 rdr. The other η-related effect is a nonaxisymmetric one: the gas cools the plasma locally and therefore suppresses the current locally, which should lead to the growth of magnetic islands with their O-point at the gas deposition location. This mechanism is involved as well in radiation driven islands, which have been invoked recently to explain density limit disruptions [START_REF] Gates | [END_REF]. To a certain extent, numerical experiments allow discriminating between the three above-mentioned mechanisms. Figure 6 shows the evolution of the magnetic energy contained in the n = 1 toroidal harmonic for three simulations. In the first one (blue), the standard η = η(T ) dependence was used. In the second one (red), we used instead η = η(T n=0 ), where T n=0 is the axisymmetric component of T , in order to suppress the non-axisymmetric ηrelated mechanism. Finally, in the third simulation (black), the T dependence of η was removed altogether by keeping the initial η profile, in order to remove both η-related effects. All three simulations are overimposed until about 0.5 ms, indicating that the early evolution is not η-related and therefore likely connected to the 3D equilibrium mechanism. However, in the third (black) simulation the n = 1 magnetic energy quickly saturates and the 2/1 island remains very small. The leading mechanisms for the 2/1 island growth are therefore ηrelated, even if the island grows from a non-η-related seed. Then, it appears that between 0.5 and 2.5 ms, the η = η(T ) case gains about one order of magnitude compared to the η = η(T n=0 ) case, showing that the non-axisymmetric η-related effect is important. This is consistent with the already mentioned experimental observation that the island O-point coincides with the gas deposition region for Ohmic plasmas [5]. Note that for neutral beam heated plasmas, the alignment is observed to degrade [5], which may be due to plasma rotation via two effects: a smearing of the nonaxisymmetric η-related effect and a drag of the island into the plasma flow.

Physics of the thermal quench

Triggering of the thermal quench

The next question is: what mechanisms trigger the TQ? This question may be addressed by analysing Figures 7 and8. Fig. 7 shows the time evolution of the magnetic energy in the different toroidal harmonics (note the logarithmic vertical scale). Fig. 8 contains poloidal cross-sections of the cosine component of ψ for the n = 1 (left), n = 2 (middle) and n = 3 (right) harmonics at different times (different rows). Both figures are for the same simulation as in Section 5. It can be seen in Fig. 7 that during the first 5 ms, magnetic perturbations are strongly dominated by the n = 1 mode. The top left plot of Fig. 8 shows that the mode has a 2/1 structure, not surprisingly. Then, at about 5 ms, a clear increase in the growth rate of the n = 2 magnetic energy is visible. This is associated to the growth of a 3/2 tearing mode, as can be seen by comparing the first and second rows, middle column, plots in Fig. 8. Around 5.8 ms, magnetic energies in higher n harmonics grow sharply. For example, a 4/3 mode can be seen to grow when comparing the second and third rows, last column, plots in Fig. 8. This leads to the TQ.

Mechanisms for the growth of the 2/1 mode have been discussed in Section 6. What about the other modes? Looking at Poincaré plots, it appears that the 3/2 mode starts growing when the inner side of the 2/1 island reaches the q = 3/2 surface. Observing the j φ,n=0 profile at different times, shown in Figure 9, it seems likely that this is due to a current profile effect. Indeed, the 2/1 island causes a local j φ,n=0 flattening, which steepens the profile radially inward (compare the black, blue and red profiles). The 3/2 mode grows when the large |j φ,n=0 | region runs across the q = 3/2 surface (i.e. between the blue and red profiles) due to the δW term in the energy principle already mentioned in Section 6. The same process then seems to take place Progress in understanding disruptions triggered by massive gas injection via 3D non-linear MHD modelling with JOREK6 with the 3/2 tearing mode, which locally flattens j φ,n=0 (red vs. magenta profiles), moving the large |j φ,n=0 | region across q = 4/3 and destabilizing the 4/3 mode. Therefore, it seems that the TQ is triggered by a kind of current profile avalanche. Note that this picture for the TQ triggering has already been proposed in [6] and [START_REF] Biskamp | Nonlinear Magnetohydrodynamics[END_REF]. The magnetic stochasticity created by the tearing modes in (roughly) the outer half of the plasma flattens the temperature profile there. In this simulation, flux surfaces however pertain inside midradius so that the central temperature cannot collapse purely from parallel thermal conduction. Convective core mixing indeed seems to play an important role in this simulation. However, as we shall see below, it is likely that tearing modes are too weak in this simulation. In other simulations where the tearing modes were excited more strongly (due to different MGI settings for example), magnetic stochasticity extended across the whole plasma. Hence, convective core mixing and the 1/1 internal kink mode may not play a large role in reality, but this is not clear at this stage.

The plasma current spike

The I p spike is a characteristic feature of disruptions. A classic explanation of its origin is that the TQ releases magnetic energy l i I 2 p (where l i is the internal inductance of the plasma) while the flux ψ at the edge of the plasma does not have the time to change significantly because the TQ is much shorter than the wall penetration time, and hence L p I p should remain constant, where L p µ 0 R[ln(8R/a) -2 + l i /2] is the self-inductance of the plasma. The consequence is that l i has to decrease and I p increase, hence the I p spike. JOREK results are well in line with this explanation, as seen in Figure 10, which displays the evolution of I p , l i , the magnetic energy inside the plasma (E mag and its approximation l i I 2 p ) and L p I p (for a different

Progress in understanding disruptions triggered by massive gas injection via 3D non-linear MHD modelling with JOREK8 simulation from above because the latter has a smaller I p spike). The fact that simulations display an I p spike in spite of having no sources answers a criticism made by Zakharov to NIMROD simulations [10] [11]. The simulated I p spike is however much smaller than the experimental one (compare Figs. 10 and1). Following the above argument, this should mean that not enough magnetic energy is released during the TQ and that l i does not decrease enough, i.e. that the current profile avalanche is not sufficiently pronounced. The strength of the avalanche may be expected to have a positive dependence on the ratio between the resistive diffusion time and the growth time of the 2/1 mode: indeed, the larger this ratio, the larger the current gradient inside the 2/1 island, and the stronger the excitation of inner modes. A stronger local cooling near q = 2 should go in this direction. The above simulation may therefore have a too weak cooling. This is also consistent with the fact that simulations over-predict the pre-TQ duration when appproaching a realistic η value. It is plausible that background impurities such as Argon remaining from previous pulses (where Ar MGI was used) [12] provide this extra cooling. The need for background impurities to match experimental observations had already been identified for He MGI simulations in Alcator C-Mod and DIII-D with NIMROD [7].

Conclusion

JOREK simulations shed light on the physics of the pre-TQ and TQ phases of a D 2 MGI-triggered disruption in an Ohmic JET plasma. The gas destabilizes (essentially via increasing η by cooling) a 2/1 tearing mode which flattens the current profile, destabilizing a 3/2 mode, which in turn destabilizes higher n modes. The energy is lost by parallel conduction along stochastic field lines and possibly also by a convective mixing of the core, but this will have to be clarified in the future. Indeed, the I p spike in the simulations is too weak, suggesting that tearing modes are not large enough, and simulations with larger tearing modes (which display a larger I p spike) produce a full stochastization of the magnetic field so that the convective core mixing plays a less important role. An interesting question which could be the object of future work is whether the TQ triggering picture found here may explain the experimental finding that the TQ is triggered at a distinct locked mode amplitude [13]. The work with JOREK on D 2 MGI simulations shall in any case be pursued, aiming for quantitative validation. For this purpose, the effect of background impurities (as well as toroidal rotation, diamagnetic effects, etc.) shall be investigated and efforts made to approach realistic parameter values, in particular for η. D 2 MGI ASDEX Upgrade simulations have been begun, for which it is easier to run at realistic parameter values. In parallel to this work, a model for MGI of other gases than D 2 has been implemented in JOREK and shall be utilized in the near future.

  Figure 1.

Figure 2 .

 2 Figure 2. Time evolution of the measured (circles) and simulated (lines -plain = ∆φMGI = 2; dashed = ∆φMGI = 0.5; dash-dotted = ∆φMGI = 100) line integrated density for the three interferometry lines of sight

Figure 3 .

 3 Figure 3. Schematic machine view from the top showing the DMV2 and interferometer chords

Figure 4 .

 4 Figure 4. Poloidal cross-sections of ne in the plane of the interferometer as simulated with JOREK for ∆φMGI = 0.5 (left) and ∆φMGI = 100 (right). Interferometry chords are also shown, with the same color code as in Fig. 2.

Figure 5 .

 5 Figure 5. Poloidal cross sections at the toroidal position of the MGI of (from left to right) Te, ne, j φ and Poincaré plots at times (from top to bottom) t = 0, t = 4.1 ms, t = 5.7 ms and t = 6.2 ms

Figure 6 .Figure 7 .

 67 Figure 6. Time evolution of the magnetic energy contained in the n = 1 component for three simulations with different models for the resistivity (see text)

Figure 8 .

 8 Figure 8. Poloidal cross-sections at the toroidal position of the MGI of (from left to right) the n = 1, n = 2 and n = 3 cosine component of ψ at times (from top to bottom) t = 4.1 ms, t = 5.1 ms and t = 5.7 ms (same simulation as in Figs. 5, 7 and 9). The color scale is the same for all plots (note the saturation for the n = 1 mode which has a large amplitude compared to the other modes).

Figure 9 .

 9 Figure 9. Toroidal current density profiles at the midplane (low field side) at different times. Note that the red and magenta profiles correspond to the last two rows of Figure 8.

Figure 10 .

 10 Figure 10. Time traces (from a different simulation than Figs. 5, 7, 9 and 8, in which the Ip spike is smaller) showing an Ip spike at the TQ (blue curve) with LpIp remaining approximately constant (black curve) while the current profile flattens (red curve) and magnetic energy is dissipated (magenta curves).

Table 1 .

 1 List of typical values of JOREK input parameters used in the present paper as well as in[3]

	Parameter	Symbol	Value in present paper	Value from [3]	Units	Comment
	Initial central resistivity	η0	3.5 × 10 -8 -3.5 × 10 -7	Same	SI (Ω • m)	Spitzer value: 2 × 10 -8
	Hyper-resistivity	η h	1 × 10 -12	5 × 10 -10	JOREK
	Initial central // thermal conductivity	κ 0	3 × 10 29	7 × 10 28	SI (/m/s)	Spitzer-Härm value: 3 × 10 29
	⊥ thermal conductivity	κ ⊥	4 × 10 19	Same	SI (/m/s)	Corresponds to χ ⊥ 1 m 2 /s
	Hyper-⊥ thermal conductivity	κ ⊥h	10 -12	Same	JOREK
	Dynamic viscosity	µ	1.4 × 10 -7	2.8 × 10 -7	SI (kg/m/s)	Corresponds to ν 1 m 2 /s
	Hyper-viscosity	µ h	1 × 10 -12	4 × 10 -10	JOREK
	// viscosity	µ	2.8 × 10 -6	2.8 × 10 -5	SI (kg/m/s)	Corresponds to ν 20 m 2 /s
	Hyper-// viscosity	µ h	10 -9	Same	JOREK
	Ion/electron particle ⊥ diffusivity	D ⊥	1.4	28	SI (m 2 /s)
	Ion/electron particle // diffusivity	D	0	2.8 × 10 4	SI (m 2 /s)
	Neutral particle diffusivity	Dn	60 -6000	2.8 × 10 4	SI (m 2 /s)

Acknowledgements

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This work was granted access to the HPC resources of Aix-Marseille University financed by the project Equip@Meso (ANR-10-EQPX-29-01). A part of this work was carried out using the HELIOS supercomputer system (IFERC-CSC), Aomori, Japan, under the Broader Approach collaboration, implemented by Fusion for Energy and JAEA, and using the CURIE supercomputer, operated into the TGCC by CEA, France, in the framework of GENCI and PRACE projects.


Van-Tan Vu
email: van-tan.vu@gipsa-lab.grenoble-inp.fr

Olivier Sename
email: olivier.sename@gipsa-lab.grenoble-inp.fr

Luc Dugard
email: luc.dugard@gipsa-lab.grenoble-inp.fr

Peter Gaspar
email: gaspar@sztaki.mta.hu

Multi objective H ∞ active anti-roll bar control for heavy vehicles

Keywords: Active anti-roll bar, H ∞ control, Multi-criteria optimization, Genetic algorithms

In the active anti-roll bar control on heavy vehicles, roll stability and energy consumption of actuators are two essential but conflicting performance objectives. In a previous work, the authors proposed an integrated model, including four electronic servo-valve hydraulic actuators in a linear yawroll model on a single unit heavy vehicle. This paper aims to design an H ∞ active anti-roll bar control and solves a Multi-Criteria Optimization (MCO) problem by using Genetic Algorithms (GAs) to select the weighting functions for the H ∞ synthesis. Thanks to GAs, the roll stability and the energy consumption are handled using a single high level parameter and illustrated via the Pareto optimality. Simulation results in frequency and time domains emphasize the efficiency of the use of the GAs method for a MCO problem in H ∞ active anti-roll bar control on heavy vehicles.

INTRODUCTION

The aim of rollover prevention is to provide the vehicle with the ability to resist overturning moments generated during vehicle maneuvers. Roll stability is determined by the height of the center of mass, the track width and the kinematic properties of the suspension. The primary overturning moment arises from the lateral acceleration acting on the center of gravity of the vehicle. More destabilizing moment arises during the cornering manoeuver when the center of gravity of the vehicle shifts laterally. The roll stability of the vehicle can be guaranteed if the sum of the destabilizing moment is compensated during a lateral manoeuver. Several schemes concerned with the possible active intervention into vehicle dynamics have been proposed. These approaches employ active anti-roll bars, active steering, active braking, active suspensions, or a combination of them [START_REF] Gaspar | The design of a combined control structure to prevent the rollover of heavy vehicles[END_REF]. The active anti-roll bar system is the most common method used to improve the roll stability of heavy vehicles. Several control design problems for active anti-roll bar systems have been investigated with many different approaches during the last decades. In (Gaspar et al., 2005a) the authors presented Linear Parameter Varying (LPV) techniques to control the active anti-roll bars, combined with an active brake control on single unit heavy vehicles. The forward velocity is considered as the varying parameter. Other works concerning the yaw-roll model on single unit heavy vehicles have dealt with optimal control (Sampson and Cebon, 2003a), robust control [START_REF] Vu | H ∞ active anti-roll bar control to prevent rollover of heavy vehicles: a robustness analysis[END_REF], and neural network control [START_REF] Boada | Active roll control using reinforcement learning for a single unit heavy vehicle[END_REF]. The H ∞ control design approach is an efficient tool for improving the performance of a closed-loop system in pre-defined frequency ranges. The key step of the H ∞ control design is the selection of weighting functions which depends on the engineer skill and experience [START_REF] Skogestad | Multivariable Feedback Control[END_REF]. In many real applications, the difficulty in choosing the weighting functions still increases significantly because the performance specification is not accurately defined i.e., it is simply to achieve the best possible performance (optimal design) or to achieve an optimally joint improvement of more than one objective (multiobjectives design). So the optimization of weighting functions to satisfy all the desired performances is still an interesting problem. In [START_REF] Hu | Systematic H ∞ weighting function selection and its application to the real-time control of a vertical take-off aircraft[END_REF] the authors proposed to consider each system, no matter how complex it is, as a combination of sub-systems of the first and second order, for which it is easy to find the good weighting functions to be used in the H ∞ control methodology. However, there is no explicit method to find these functions in the general case. The usual way is to proceed by trial-and-error. Recently, the idea to use an optimization tool was proposed in [START_REF] Alfaro-Cid | Optimisation of the weighting functions of an H ∞ controller using genetic algorithms and structured genetic algorithms[END_REF]. The choice of GAs seems natural because their formulation is well suited for this type of problems [START_REF] Do | Multi-objective optimization by genetic algorithms in H ∞ /LPV control of semi-active suspension[END_REF].

Based on the integrated model presented in (Vu et al., 2016a), this paper proposes an H ∞ control for active anti-roll bar, the GAs method is used to solve the Multi-Criteria Optimization (MCO) problem for the H ∞ synthesis. The latter work is here extended providing two new main contributions:

• We design here an H ∞ controller for active anti-roll bar system on the integrated model for the single unit heavy vehicle. The aim is to improve the roll stability of the heavy vehicle. The normalized load transfers and the limitation of the input currents generated by the controllers are considered. • The GAs method is applied to define the weighting functions of the H ∞ controller for the active anti-roll bar system on the integrated model. Thanks to GAs, the conflicting objectives between the normalized load transfers and the input currents are handled using only one single high level parameter. This paper is organised as follows: Section 2 gives the integrated model for a single unit heavy vehicle. Section 3 presents the MCO problem of active anti-roll bar control. Section 4 illustrates the H ∞ robust control synthesis to improve roll stability of heavy vehicles. In section 5, GAs method is used for MCO in the H ∞ anti-roll bar control. Section 6 shows the simulation results in the frequency and time domains. Finally, some conclusions are drawn in section 7.

INTEGRATED MODEL FOR HEAVY VEHICLES

The proposed integrated model includes four Electronic Servo-Valve Hydraulic (ESHV) actuators (two at the front axle and two at the rear axle) in a linear single unit heavy vehicle yawroll model [START_REF] Gaspar | The design of an integrated control system in heavy vehicles based on an lpv method[END_REF]. The control signal is the electrical current u opening the electronic servo-valve, the output is the force F act generated by the hydraulic actuator. The symbols of the integrated model are found in Table 1. In the linear single unit heavy vehicle yaw-roll model, the differential equations of motion, i.e., the lateral dynamics, the yaw moment, the roll moment of the sprung mass, the roll moment of the front and the rear unsprung masses, are formalized in the equations (1):

                                             mv( β + ψ) -m s h φ = F y f + F yr -I xz φ + I zz ψ = F y f l f -F yr l r (I xx + m s h 2 ) φ -I xz ψ = m s ghφ + m s vh( β + ψ) -k f (φ -φ t f ) -b f ( φ -φt f ) + M AR f + T f -k r (φ -φ tr ) -b r ( φ -φtr ) + M ARr + T r -rF y f = m u f v(r -h u f )( β + ψ) + m u f gh u f .φ t f -k t f φ t f +k f (φ -φ t f ) + b f ( φ -φt f ) + M AR f + T f -rF yr = m ur v(r -h ur )( β + ψ) -m ur gh ur φ tr -k tr φ tr +k r (φ -φ tr ) + b r ( φ -φtr ) + M ARr + T r (1)
In (1) the lateral tyre forces F y;i in the direction of velocity at the wheel ground contact points are modelled by a linear stiffness as:

F y f = µC f α f F yr = µC r α r (2)
with tyre side slip angles:

             α f = -β + δ f - l f ψ v α r = -β + l r ψ v (3)
The moment of passive anti-roll bar impacts the unsprung and sprung masses at the front and rear axles as follows:

               M AR f = 4k AO f t A t B c 2 φ -4k AO f t 2 A c 2 φ u f M ARr = 4k AOr t A t B c 2 φ -4k AOr t 2 A c 2 φ ur (4)
where k AO f , k AOr are respectively the torsional stiffnesses of the anti-roll bar at the front and rear axles, t A is half the distance of the two suspensions, t B is half the distance of the chassis, c is the length of the anti-roll bars' arm.

The torque generated by the active anti-roll bar system at the front axle is now determined by: T f = 2l act F act f = 2l act A p ∆ P f (5) and the torque generated by the active anti-roll bar system at the rear axle is:

T r = 2l act F actr = 2l act A p ∆ Pr (6) where ∆ P f and ∆ Pr are respectively the difference of pressure of the hydraulic actuator at the front and rear axles. The equations of these electronic servo-valve actuators are given by (7): Effective bulk modulus of the oil τ Time constant of the servo-valve

K v Servo-valve gain X v f,r
Spool valve displacements at the front/rear axle

                                         V t 4β e ∆P f + (K P + C tp )∆ P f -K x X v f +A p l act φ -A p l act φu f = 0 Ẋv f + 1 τ X v f - K v τ u f = 0 V t 4β e ∆Pr + (K P + C tp )∆ Pr -K x X v f +A p l act φ -A p l act φur = 0 Ẋvr + 1 τ X vr - K v τ u r = 0 (7)
Defining the state vector:

x = β ψ φ φ φ u f φ ur ∆ P f X v f ∆ Pr X vr T
The motion differential equations ( 1)-( 7) can be rewritten in the LTI state-space representation as:

ẋ = A.x + B 1 .w + B 2 .u z = C.x + D 1 .w + D 2 .u (8)
where A, B 1 , B 2 , C, D 1 , D 2 are the model matrices of appropriate dimensions.

The exogenous disturbance (steering angle) is:

w = δ f T
and the control inputs (input currents):

u = u f u r T

MULTI-CRITERIA OPTIMIZATION OF ACTIVE ANTI-ROLL BAR CONTROL

Multi-criteria optimization and Pareto-optimal solutions

Multi-Criteria Optimization (MCO) can be described in mathematical terms as follows [START_REF] Ehrgott | Multicriteria optimization[END_REF]:

min x∈S F(x) = f 1 (x), f 2 (x), ..., f n (x) (9)
where n > 1 and S is the set of constraints defined above.

The space in which the objective vector belongs is called the objective space, and the image of the feasible set under F is called the attained set. In the following such a set will be denoted by C = {y ∈ R n : y = f (x), x ∈ S }. The scalar concept of optimality does not apply directly in the multi-criteria setting.

Here the notion of Pareto optimality is introduced. Essentially, a vector x * ∈ S is said to be Pareto optimal for a multi-criteria problem if all other vectors x ∈ S do have a higher value for at least one of the objective functions f i , with i = 1, ..., n, or have the same value for all the objective functions.

There are many formulations to solve the problem ( 9) such as weighted min-max method, weighted global criterion method, goal programming methods... [START_REF] Marler | Survey of multi-objective optimization methods for engineering[END_REF]) and references therein. Here, one uses a particular case of the weighted sum method, where the multi-criteria functions vector F is replaced by the convex combination of objectives:

min x∈S J = n i=1 α i f i (x), n i=1 α i = 1 (10)
The vector α = (α 1 , α 2 , ..., α n ) represents the gradient of function J.

Control objective, and MCO for H ∞ active anti-roll bar control

The main objective of the active anti-roll bar control system is to maximize the roll stability of the vehicle to prevent a rollover phenomenon in an emergency. Two main criteria are commonly used to assess the roll stability of the heavy vehicle:

• The normalized load transfer R f,r at the two axles, defined as follows [START_REF] Hsun-Hsuan | Active roll control for rollover prevention of heavy articulated vehicles with multiple-rollover-index minimisation[END_REF]:

R f = ∆F z f F z f , R r = ∆F zr F zr (11)
where F z f is the total axle load at front axle and F zr at rear axle. ∆F z f and ∆F zr are respectively the lateral load transfers at the front and rear axles, which can be given by:

∆F z f = k u f φ u f l w , ∆F zr = k ur φ ur l w (12)
where k u f and k ur are the stiffness of the tyres, φ u f and φ ur are the roll angles of the unsprung masses at the front and rear axles, l w the half of the vehicle's width.

The normalized load transfer R f,r = ±1 value corresponds to the largest possible load transfers. The roll stability is achieved by limiting the normalized load transfers within the levels corresponding to wheel lift-off. • The roll angles between the sprung and unsprung masses (φ -φ u ), give the maximum stabilizing moment of the active anti-roll bar system to be increased. They should stay within the limits of the suspension travel 7 -8deg [START_REF] Sampson | Active roll control of single unit heavy road vehicles[END_REF]). As mentioned above, the objective of the active anti-roll bar control system is to improve the roll stability of heavy vehicles. However, such a performance objective must be balanced with the energy consumption of the anti-roll bar system due to the input current entering the electronic servo-valve of the actuators. Therefore the objective function is selected as follows:

f = α f Normalized-load-trans f er + (1 -α) f Control-cost (13) The vector α = [0 ÷ 1]
is the gradient of function f . When α moves to 0, the optimal problem focusses on minimizing input currents. And conversely, when α moves to 1, the optimal problem focusses on minimizing normalized load transfers. In the objective function ( 13), f Normalized-load-trans f er and f Control-cost are performance indices corresponding to the normalized load transfers and input currents at the two axles, which are defined as follows:

                         f Normalized-load-trans f er = 1 2          1 T T 0 R 2 f (t)dt + 1 T T 0 R 2 r (t)dt          f Control-cost = 1 2             1 T T 0 u 2 f (t)dt 1 T T 0 u 2 f (t) max dt + 1 T T 0 u 2 r (t)dt 1 T T 0 u 2 r (t) max dt             ( 14
)
where u f,rmax are defined when the optimal problem focusses only on the normalized load transfers (i.e., the input currents are then not considered in the optimisation problem). In that case, α = 1 and f = f Normalized-load-trans f er .

The MCO problem represented by the equation ( 13) can not be resolved directly in the synthesis of H ∞ controller. Thus, summarizing the implementation is done in this paper is described as Fig 2 . The generalized plant includes the integrated model (section 2) and the weighting functions. The controller is synthesised in term of H ∞ method (section 4). The conflicting objective between roll stability and energy consumption is the computation of the close-loop performance (MCO problem in section 3). Depending on the purpose of the MCO problem, the weighting functions are appropriately selected by GAs (section 5). The optimal parameters obtained from GAs are sent to the weighting functions to calculate the characteristics of the closed-loop system (section 6).

H ∞ ACTIVE ANTI-ROLL BAR CONTROL TO IMPROVE ROLL STABILITY OF HEAVY VEHICLES

Background on H ∞ control

The interested reader may refer to [START_REF] Skogestad | Multivariable Feedback Control[END_REF], [START_REF] Scherer | Linear matrix inequalities in control[END_REF] for detailed explanations on H ∞ control design. The aim is to design a controller K(s) that reduces the signal transmission path from disturbances d to performance outputs z and also stabilizes the closed-loop system. The H ∞ problem is to find K which minimizes γ such that F l (P, K) ∞ < γ (18) By minimizing a suitably weighted version of (18), the control aim is achieved, as presented below. Figure 4 shows the closed-loop structure of an H ∞ control designed for the active anti-roll bar system on a single unit heavy vehicle, using ESVH actuators. In the diagram, the feedback structure includes the nominal model G, the controller K, the performance output z, the control input u, the measured output y, the measurement noise n. The steering angle δ f is the disturbance signal, which is set by the driver. The weighting functions W δ , W z , W n are presented below.

H ∞ control design for active anti-roll bar system

According to Figure 4, the concatenation of the linear model ( 8) with the performance weighting functions lead to the state space representation of P(s):

        ẋ z y         =         A B 1 B 2 C 1 D 11 D 12 C 2 D 21 D 22                 x w u         (19)
with the exogenous input:

w = δ f n
the control input:

u = u f u r T
where u f and u r are the input current at the two axles. the performance output vector:

z = u f u r R f R r a y
T and the measured output vector:

y = a y φ T
The input scaling weight W δ normalizes the steering angle to the maximum expected command. It is selected as W δ = π/180, which corresponds to a 1 0 steering angle command. The weighting function W n is selected as a diagonal matrix, which accounts for small sensor noise models in the control design. The noise weights are chosen as 0.01(m/s 2 ) for the lateral acceleration and 0.01( 0 /sec) for the derivative of roll angle φ [START_REF] Gaspar | The design of a combined control structure to prevent the rollover of heavy vehicles[END_REF].

The weighting functions matrix W z represents the performance output, W z = diag[W zu , W zR , W za ]. The purpose of the weighting functions is to keep the control inputs, normalized load transfers and lateral acceleration as small as possible over the desired frequency range. These weighting functions can be considered as penalty functions, that is, weights should be large in the frequency range where small signals are desired and small where larger performance outputs can be tolerated. The weighting function W zu is chosen as

W zu = diag[W zu f , W zur ],
corresponding to the input currents at the front and rear axles, and are chosen as:

W zu f = 1 Z 1 ; W zur = 1 Z 2 (20)
The weighting function W zR is chosen as

W zR = diag[W zR f , W zRr ],
corresponding to the normalized load transfers at front and rear axles, and are selected as:

W zR f = 1 Z 3 ; W zRr = 1 Z 4 (21)
The weighting function W za is selected as:

W za = Z 51 Z 52 s + Z 53 Z 54 s + Z 55 (22)
Here, the weighting function W za corresponds to a design that avoids the rollover situation with the bandwidth of the driver in the frequency range up to more than 4rad/s. This weighting function will directly minimize the lateral acceleration when it reaches the critical value, to avoid the rollover.

As said before, the key step of the H ∞ control design is how to select the weighting function. This is not an easy task even with good engineers, because it depends on the engineer skill and experience. From equations ( 20) -( 22), Z i and Z 5, j are constant parameters. The following variables are to be defined: Z 1 , Z 2 , Z 3 , Z 4 , Z 51 , Z 52 , Z 53 , Z 54 , Z 55 . In the next section, the GAs method will be used to find these variables, suited for the MCO problem.

USING GENETIC ALGORITHMS FOR MULTI-CRITERIA OPTIMIZATION IN H ∞ ANTI-ROLL BAR CONTROL

This section introduces the MCO problem for the H ∞ active anti-roll bar control on heavy vehicles, which is solved by using the GAs method.

Genetic Algorithms

A Genetic Algorithms developed by J.H. Holland [START_REF] Holland | Adaptation in natural and artificial systems: an introductory analysis with applications to biology, control, and artificial intelligence[END_REF] is a model of machine learning, which derives its behavior from a metaphor of the processes of evolution in nature. GAs are executed iteratively on a set of coded chromosomes, called a population, with three basic genetic operation: selection, crossover and mutation. Each member of the population, called a chromosome (or individual) is represented by a string. GAs use only the objective function information and probabilistic transition rules for genetic operations. The primary operation of GAs is the crossover. The crossover happens with a probability of 0.9 and the mutation happens with a small probability 0.095. The basic flowchart of GAs is shown in Fig 5.

Solving multi-criteria optimization by genetic algorithms

From the objective function in ( 13), the MCO problem for the H ∞ active anti-roll bar control can be defined as:

min p∈P F(p), F(p) := f Normalized-load-trans f er , f Control-cost T P := p = [Z 1 , Z 2 , Z 3 , Z 4 , Z 51 , Z 52 , Z 53 , Z 54 , Z 55 ] T ∈ R | p l ≤ p ≤ p u (23)
where F(p) is the vector of objectives, p denotes the vector of the weighting function parameters, and p l , p u represent the lower and upper bounds of the weighting function selection.

The lower and upper bounds of the weighting function parameters are given in Table 2. Besides the minimization of the objective function from equations ( 13) and ( 23), we also have to account for the limitations of the normalized load transfers, roll angle of suspensions, as well as the input currents at each axle. These limitations are considered as the optimal conditions (binding conditions) shown in the Table 3. The proposed weighting function optimization procedure for the H ∞ active anti-roll bar control synthesis is as follows:

| φ -φ u f | < 7 deg 2 | φ -φ ur | < 7 deg 3 | R f | < 1 - 4 | R r | < 1 - 5 | u f | < 20 mA 6 | u r | < 20 mA
•

Step 1: Initialize with the weighting functions (it depends on the engineer skill and experience), the vector of weighting functions selected as p = p 0 . • Step 2: Select lower bound, upper bound, scaling factor, offset and start point. • Step 3: Format optimal algorithm, select the vector of objectives with the variation of switch value from 0 to 1 and then solve the minimization problem. 6. SIMULATION RESULTS

Optimization results

Thanks to the GAs method, Table 4 gives a synthesis of the values of the variables Z i , Z 5 j in six cases for α = [0; 0.25; 0.5; 0.7; 0.9; 1], as explained in (13). When α = 0, f = f Control-cost , the optimal problem focusses only on the input currents and when α = 1, f = f Normalized-load-trans f er , the optimal problem focusses only on the normalized load transfers. To evaluate the optimization procedure, the simulation results in both frequency and time domains are done and compared for five cases: passive anti-roll bar (Vu et al., 2016a) and H ∞ active anti-roll bar with α = [0; 0.5; 0.9; 1] .

Evaluation of optimization results in frequency domain

In this section, the frequency response of the heavy vehicle is shown in the nominal parameters case of the single unit heavy vehicle being considered, characterized by the forward velocity V at 70Km/h and the road adhesion coefficient µ = 1 (see, [START_REF] Gaspar | The design of a combined control structure to prevent the rollover of heavy vehicles[END_REF] and Vu et al. (2016a)). Figures 7 and8 show the transfer function magnitude of the normalized load transfers at the two axles To assess the roll stability of the heavy vehicle using the four H ∞ active anti-roll bar controllers. The reduction of the magnitude of transfer functions compared with the passive antiroll bar case is considered at 10 -2 rad/s and at 2rad/s as:

R f,r δ f .
λ (X) = X active δ f - X passive δ f ( 24
)
where X is the variable of interest, which are the normalized load transfers R f,r .

Figure 9 shows the reduction of the magnitude of transfer functions of the normalized load transfer compared with the passive anti-roll bar case at 10 -2 rad/s and at 2rad/s. We can see that at 10 -2 rad/s the controller with α = 0 decreases the roll stability, meanwhile, when the value of α increases, the roll stability of the heavy vehicle increases. The curves are very different. At 2rad/s the transfer functions are not so different. This will be explained in further studies. δ f : when α increases (the MCO problem focusses on the minimizing the normalized load transfers), the total input currents also increase. It's proven for the normalized load transfer and the input current that they are two conflicting performance objectives. Thus through the MCO problem, we can choose the weighting functions to enhance the roll stability of the heavy vehicle in the low frequency range as well as in the high frequency range up to over 4rad/s which is the limited bandwidth of the driver [START_REF] Gaspar | The design of a combined control structure to prevent the rollover of heavy vehicles[END_REF].

Evaluation of optimization results in time domain

In this section, the cornering responses of a single unit vehicle can be seen. The steering angle (disturbance) applied in the simulation is a step signal as in Fig 12. [START_REF] Gaspar | The design of a combined control structure to prevent the rollover of heavy vehicles[END_REF].

The forward velocity of the heavy vehicle continuously varies during operation, especially in the case of an emergency. The rollover of the heavy vehicle often occurs with a forward velocity within 60 to 110Km/h. In Figures 13 -14 we consider the forward velocity of the heavy vehicle up to 160Km/h in order to evaluate the roll stability when the normalized load transfers reach its limitations.

From Figure 13, we can see that the maximum absolute value of normalized load transfers at the front axle reaches the limit "1" in the case of α = [0; 0.5; 0.9; 1] where the forward velocities are respectively 74, 110, 139, 144Km/h. Note that in the case of the passive anti-roll bar, we get 92Km/h for the forward velocity. Considering Figure 14, the maximum absolute value of the normalized load transfers at the rear axle reaches the limit "1 in the case of α = [0; 0.5; 0.9; 1] where the forward velocities are respectively 66, 84, 129, 139Km/h. Note that in the case of the passive anti-roll bar, we get 80Km/h for the forward velocity.

From Figures above, we can conclude that the risk of the rollover of heavy vehicles is reduced thanks to the multi objective H ∞ controller.

CONCLUSION

In this paper, the integrated model of a single unit heavy vehicle including four ESVH actuators is used to develop a linear H ∞ control scheme maximizing its roll stability in order to prevent rollover. The normalized load transfers and the limitations of the input currents are considered in the development.

A weighting function optimization procedure using GAs for H ∞ active anti-roll bar control on the integrated model has also been proposed. The conflicting objectives between the normalized load transfers and input currents are handled using only one high level parameter, which is a great advantage to solve the multi-objective control problem. The simulation results in frequency and time domains have shown the efficiency of GAs in finding a suitable controller to satisfy the MCO problem. Even if a LTI controller seems to performs reasonably well here, the comparison with an LPV controller (scheduled by the vehicle velocity) will be of interest for future works.

  By using various sets of α, one can generate several points in the Pareto set. The shape of the Pareto surface indicates the nature of the tradeoff between the different objective functions. An example of a Pareto curve is reported in Fig 1, where all the points between ( f 2 ( x), f 1 ( x)) and ( f 2 ( x), f 1 ( x)) define the Pareto front. These points are called non-inferior points.

Fig. 1 .

 1 Fig. 1. Example of a Pareto curve.

Fig. 2 .

 2 Fig. 2. Controller optimization for H ∞ active anti-roll bar using Genetic Algorithms. The H ∞ control problem is formulated according to the generalized control structure shown in Fig 3.

Fig

  Fig. 3. Generalized control structure. with P partitioned as z y = P 11 (s) P 12 (s) P 21 (s) P 22 (s) d u (15) and u = K(s).y (16) which yields z d = F l (P, K) := [P 11 + P 12 K[I -P 22 K] -1 P 21 ] (17)

Fig. 4 .

 4 Fig. 4. Closed-loop structure of an H ∞ active anti-roll bar control.

Fig. 5 .

 5 Fig. 5. Basic flowchart of Genetic Algorithm.

•

  Step 4: Select the individuals, apply crossover and mutation to generate a new generation: p = p new . • Step 5: Evaluate the new generation by comparing with the binding conditions. If the criteria of interest are not satisfied, go to step 3 with p = p new ; else, stop and save the best individual: p opt = p new .

Fig. 6 .

 6 Fig. 6. The Pareto frontier for the active anti-roll bar on heavy vehicles using ESVH actuators.

Fig. 8 .

 8 Fig. 7. Transfer function magnitude of the normalized load transfer at the front axle R f δ f .

Fig. 9 .

 9 Fig. 9. Reduction of the magnitude of transfer functions of the normalized load transfers at the two axles compared with the passive anti-roll bar case (see (24)). Figures 10 and 11 show the transfer function magnitude of the input currents at the two axles u f,r

  Fig. 10. Transfer function magnitude of the input current at the front axle u f δ f .

Fig. 11 .

 11 Fig. 11. Transfer function magnitude of the input current at the rear axle u r δ f .

Fig. 12 .

 12 Fig. 12. Steering angle in the cornering manoeuvre[START_REF] Gaspar | The design of a combined control structure to prevent the rollover of heavy vehicles[END_REF]).

Fig. 13 .

 13 Fig. 13. Effect of the forward velocity on the normalized load transfer: front axle R f .

Fig. 14 .

 14 Fig. 14. Effect of the forward velocity on the normalized load transfer: rear axle R r .

Table 1 .

 1 Symbols of the integrated model (seeGaspar et al. (2005a),Vu et al. (2016a)).

	Symbols Description
	m s	Sprung mass
	m u f,r	Unsprung mass on the front/rear axle
	m	The total vehicle mass
	v	Forward velocity
	h	Height of sprung mass from roll axis
	h u,i	Height of unsprung mass from ground
	r	Height of roll axis from ground
	a y	Lateral acceleration
	β	Side-slip angle at center of mass
	ψ	Heading angle
	ψ	Yaw rate
	α	Side slip angle
	φ	Sprung mass roll angle
	φ u f,r	Unsprung mass roll angle at the front/rear axle
	δ f	Steering angle
	u f,r	Control current at the front/rear axle
	C f,r	Tire cornering stiffness on the front/rear axle
	k f,r	Suspension roll stiffness on the front/rear axle
	b f,r	Suspension roll damping on the front/rear axle
	k t f,r	Tire roll stiffness on the front/rear axle
	I xx	Roll moment of inertia of sprung mass
	I xz	Yaw-roll inertial of sprung mass
	I zz	Yaw moment of inertia of sprung mass
	l f,r	Length of the front axle/rear axle from the CG
	l w	Half of the vehicle width
	µ	Road adhesion coefficient
	A P	Area of the piston
	K x	Valve flow gain coefficient
	K P	Total flow pressure coefficient
	C tp	Total leakage coefficient
	V t	Total volume of trapped oil
	β e

Table 2 .

 2 Lower and upper bounds of the weighting functions.

	Bounds	Wzu f Z1	Wzur Z2	WzR f Z3	WzRr Z4	Z 51	Z 52	Wza Z 53	Z 54	Z 55
	Lower bound	0.001	0.001	0.1	0.1	0.5	1 3000	1	10	0.001
	Upper bound	10	10	100	100	100	10	900 1000	20

Table 3 .

 3 Binding conditions.

	No	Note	Maximum value Unit
	1

Table 4 .

 4 Optimization results for the weighting functions of H ∞ active anti-roll bar.

	Controllers	W zu f Z 1	W zur Z 2	W zR f Z 3	W zRr Z 4	Z 51	Z 52	W za Z 53	Z 54	Z 55
	α = 0	0.060	0.020	0.100	0.965	0.673	0.948	1.063	972.212	0.855
	α = 0.25	0.057	0.052	0.51	0.863	0.863	0.664	155.627	651.707	0.573
	α = 0.5	0.099	0.0773	1.403	0.217	0.812	0.813	88.666	407.658	1.001
	α = 0.7	0.057	0.066	0.412	0.221	0.832	0.514	139.609	357.401	1.901
	α = 0.9	0.066	0.072	0.616	0.482	0.724	0.492	202.316	455.747	0.544
	α = 1	0.07	0.090	0.655	0.305	0.545	0.245	444.397	839.299	0.163