A Jaafar
email: jaafar@laplace.univ-tlse.fr

C Turpin
email: turpin@laplace.univ-tlse.fr

X Roboam
email: roboam@laplace.univ-tlse.fr

E Bru
email: bru@laplace.univ-tlse.fr

O Rallieres
email: rallieres@laplace.univ-tlse.fr

Energy management of a hybrid system based on a fuel cell and a Lithium Ion battery: experimental tests and integrated optimal design

Keywords: Fuel Cell, Battery, Lithium Ion, Energy Management Strategy, Optimal Design Nomenclature FC: Fuel Cell BAT: Battery SFEM: Symmetric Frequency Energy Management AFEM: Asymmetric Frequency Energy I_FC: fuel cell current I_load: load current Ich_BAT

The optimal design of multisource systems, hybrid systems in particular, requires an adequate choice of the energy management strategy. This latter usually impacts source sizing and lifetime. The present paper deals with an energy management approach based on a frequency sharing of the mission. Firstly, the limits of a symmetric frequency energy management are presented in the case of a hybrid system associating a fuel cell with a Li-Ion battery. Subsequently, an original energy "asymmetric" management strategy for the optimal sizing of this association is presented. This strategy is then tested on the "Hydrogen" platform at the LAPLACE research laboratory. Finally, the two energy management strategies are compared in the context of an integrated design by optimization; the asymmetric strategy offers significant gains in terms of system weight, which is important for embedded applications.

Introduction

Considering major energy challenges of this new century and in the context of fossil fuel depletion and greenhouse gas emissions, hydrogen is a promising solution for the storage of renewable energies and can be seen as "a potential energy vector for the future" [START_REF] Péra | Electrochemical components[END_REF]. Its transformation into electrical energy is ensured by a Fuel Cell (FC). Unlike storage components (batteries, ultra-capacitors, flywheels) where energy and power are closely linked in the same component, FCs offer energy / power decoupling. FC system can then be considered as an energy source (with a high specific energy range with respect to the batteries) whose autonomy is related to the hydrogen tank size. As a consequence, the delivered power only depends on FC stack sizing.

Therefore, fuel cells are utilised in a wide range of applications, particularly in transport or embedded systems [START_REF] Chan | The State of the Art of Electric, Hybrid, and Fuel Cell Vehicles[END_REF], [START_REF] Cropper | Fuel Cells: A survey of current developments[END_REF] where weight and volume constraints are the key drivers, but also in stationary field [START_REF] Hamelin | Dynamic behavior of a PEM fuel cell stack for stationary applications[END_REF] for which cost and efficiency are more concerned. In all these fields of application, the lifetime of such devices is also questioned and may be greater than 10,000 hours for constant operation. However, as soon as the electrical power draw becomes intermittent, the FC lifetime is significantly reduced. To overcome this difficulty, major applications rely on the hybridization of sources associating the FC with storage devices, while also defining an energy management strategy which optimizes power sharing between the FC and storage devices. In [START_REF] Büchi | Electrical Drivetrain Concept with Fuel Cell System and Supercapacitor Results of the "HYLIGHT®"-vehicle[END_REF], [START_REF] Martínez | Practical Control Structure and Energy Management of a Test bed Hybrid Electric Vehicle[END_REF], [START_REF] Drolia | An approach to connect ultracapacitor to fuel cell powered electric vehicle and emulating fuel cell electrical characteristic, using switched mode converter[END_REF] and [START_REF] Davat | Fuel cell-based hybrid systems[END_REF], authors have studied FC based hybrid structures with management strategies for electric vehicle applications; a double hybridization of FC source with ultra-capacitors and batteries is proposed in [START_REF] Thounthong | Energy management of fuel cell/battery/supercapacitor hybrid power source for vehicle applications[END_REF]. FC powered hybrid systems have also been studied for railway traction applications [START_REF] Miller | System design of a large fuel cell hybrid locomotive[END_REF] or in more electric aircraft [START_REF] Arregui | Theoretical study of a power generation unit based on the hybridization of a fuel cell stack and ultracapacitors[END_REF]; in this field of application, system weight has to be minimized as proposed in this paper. In [START_REF] Morin | Direct hybridization of Fuel cell -Ultracapacitors[END_REF], authors propose an original direct association of fuel cell devices with ultra-capacitor storage without any power electronic interface; the issue is here to enhance the life duration of fuel cells by means of this direct hybridization.

In order to associate efficiently several sources in a hybrid system, various types of energy management strategies can be applied. For example, "slope power sharing" [START_REF] Bang | New high power -high ratio non isolated DC-DC boost converter for Fuel cell applications[END_REF], "band power sharing" [START_REF] Büchi | Electrical Drivetrain Concept with Fuel Cell System and Supercapacitor Results of the "HYLIGHT®"-vehicle[END_REF], fuzzy logic based power sharing [START_REF] Martínez | Practical Control Structure and Energy Management of a Test bed Hybrid Electric Vehicle[END_REF], [START_REF] Papadimitriou | Integration of a hybrid fuel cell-battery system to a distribution grid[END_REF], "optimal planning" based power sharing [START_REF] Rajaei | Control Strategies for Hybrid Electric Vehicles: Evolution, Classification, Comparison, and Future Trends[END_REF], [START_REF] Zhu | A four-step method to design an energy management strategy for hybrid vehicles[END_REF], [START_REF] Feroldi | Energy Management Strategies based on efficiency map for Fuel Cell Hybrid Vehicles[END_REF] and "frequency power sharing" [START_REF] Chapoulie | Modelling of an Electric Vehicle Including Ultracapacitors with SABER[END_REF], [START_REF] Arregui | Theoretical study of a power generation unit based on the hybridization of a fuel cell stack and ultracapacitors[END_REF], [START_REF] Jaafar | Sizing and energy management of a hybrid locomotive based on flywheel and accumulators[END_REF] may be referred to.

In the present paper, we have studied and compared two energy management strategies, each of them being based on a frequency approach; this kind of strategy is applied to a particular hybrid system composed of a Fuel Cell device connected to a Li-Ion battery system.

Classically, a Symmetric Frequency Energy Management (SFEM, described in sub section 2.1) is proposed as power management strategy; its main characteristic is related to the symmetry of power profile drawn from the storage device, limiting therefore the decrease of SOC during the system operation [START_REF] Turpin | Systemic Design Methodologies for Electrical Energy Systems[END_REF]. However, contrarily to the case of ultra-capacitor storage, the limitation of the SFEM strategy is mainly due to the asymmetry of charge vs discharge current capacities for Lithium Ion storage technologies. This asymmetry may lead to an oversizing of the storage device in order to be able to accept maximum charge current constraints then increasing the system weight. Due to this drawback, the SFEM has been compared with an original management strategy called « Asymmetric Frequency Energy Management » (AFEM) which is specific to the Fuel Celllithium Ion storage association. Indeed, the asymmetry of charge/discharge current constraints on the accumulator is in accordance with the asymmetrical capabilities of fuel cell systems as described in sub section 2.2. The AFEM strategy is compared with the previous one by simulations (section 3) and has been experimentally tested on the « Hydrogen » platform of the LAPLACE lab as presented in section 4.

In the last section, an optimal design of the hybrid association integrating both (SFEM AFEM) strategies is proposed. This approach is especially dedicated to the field of embedded systems, especially for more electric aircraft for which geometrical constraints are the key issue. The system weight is then targeted as the most important design objective to be minimized together with system losses. Finally, the advantages of this original asymmetrical energy management strategy on the system weight are assessed in the context of a systemic design approach by optimization.

Energy management of the hybrid system

Symmetric frequency energy management

"Frequency" energy management is based on satisfying energy and power constraints of each source given its specific power/energy capability. It ensures the compatibility between the frequency components of the mission and the intrinsic frequency capacities of the different sources. In the case of the studied hybrid system, Symmetric Frequency Energy Management (SFEM) consists in providing the low frequency part of the mission by the FC.

The high frequency part of the mission is then devoted to the battery. The implementation of this energy management strategy is rather simple and consists in using a low pass filter as shown in Figure 1. We denote by Fg the cutoff frequency of the filter. Let us note that this management strategy based on a frequency sharing does not manage regenerative load currents. Indeed, the FC system being non power reversible, potential regenerative currents should be necessarily accepted by the battery. Consequently, the battery profile would not remain symmetric. This issue has to be taken into account for electric vehicle application. However, our privileged class of systems deals with embedded systems as for aircraft electric networks for which only positive global power consumption are considered. the mission (I_FC) and the battery provides the high frequency part (I_BAT). We can observe that according to the sign convention chosen in the study, positive currents correspond to discharge while negative currents deal with charge phases.

The battery mission profile (I_BAT) obtained in Figure 3.a shows that the charge phases have the same maximum current (in absolute value) as the discharge phases. Furthermore, based on the symmetry of charge vs discharge phases, this mission is characterized by a low energy content (simply due to the storage system losses) compared to the energy required for the global load mission (I_load). This symmetry of charge vs discharge phases has two main consequences:

-On the one hand, only losses in the batterythe chopperthe cable association have to be compensated as these latter provoke a very slow decrease of the battery SOC. A supplementary current (Ich_BAT) has to be provided by the fuel cell to compensate this slow decrease. Several strategies are possible to control the battery SOC or to limit its slow decrease as explained in [START_REF] Arregui | Theoretical study of a power generation unit based on the hybridization of a fuel cell stack and ultracapacitors[END_REF] for a fuel cellultra capacitor hybrid system. As illustrated on Figure 4, a battery recharging open loop control may be added. Note that this compensation loop has not actually been implemented in experiments based on SFEM strategy. Indeed, during the time interval of the mission profile (less than 4 hours) the very slow decrease of SOC (only due to losses) makes unnecessary to implement the SOC compensation loop. -On the other hand, considering the natural dissymmetry between the Li Ion battery charge and discharge current capacities (conversely to the case of ultra-capacitors), the sizing of the storage device is clearly constrained by the charge current. For example, the manufacturer of our Li-Ion battery recommends a maximum discharge current of 90 A and a maximum charge current of only 6 A. Consequently, the battery necessary to satisfy the I_BAT mission is oversized in terms of energy and power. This is solely due to the maximum charge currents; this issue is emphasized in section 5 where optimal design constraints are studied. To solve this problem, in the next section a proposal related to Asymmetric Frequency Energy Management strategy (AFEM) is outlined.

Asymmetric frequency energy management

Asymmetric Frequency Energy Management (AFEM) principle implies providing the low frequency part of the mission by the FC as for SFEM, but only during phases characterized by a positive gradient of current. Subsequently, the battery provides a discharge current which corresponds to the rest of the mission, i.e., the high frequency components of the mission power during positive gradients. Conversely, the battery is not charged during the negative gradient of the mission profile, this load being fully achieved by the FC system.

For the positive gradients, the major risk for the FC is the gas starvation during the transitory periods due to the time response of the gas supplies [START_REF] Borup | Scientific Aspects of Polymer Electrolyte Fuel Cell Durability and Degradation[END_REF], [START_REF] Lin | Investigation of dynamic driving cycle effect on the degradation of proton exchange membrane fuel cell by segmented cell technology[END_REF], [START_REF] Yousfi-Steiner | A review on polymer electrolyte membrane fuel cell catalyst degradation and starvation issues: Causes, consequences and diagnostic for mitigation[END_REF]. So, because of too small gas flows during these transient periods, the carbon electrodes will be a little bit consumed (CO2 production). By repeating a lot of times the positive gradients, the degradation of the electrodes is then accelerated.

For the negative gradients, the major risk could be the membrane drying due to the time response of the gas supplies again (the gas flows are too high during the transient periods) but the time constant for the drying is more than several minutes, so this drying phenomenon cannot occur in practice [START_REF] Barbir | Relationship between pressure drop and cell resistance as a diagnostic tool for PEM fuel cells[END_REF], [START_REF] Fouquet | Model based PEM fuel cell state-of-health monitoring via ac impedance measurements[END_REF]. The negative gradients are then not degrading for the FC.

According to this energy management strategy, the mission of the battery would not include charge currents. Consequently, contrarily to the SFEM case, the battery charge would not naturally be ensured. A battery slow charging loop respecting manufacturer's recommendations should be added. The battery may be charged with a current (Ich_BAT) which varies linearly according to its state of charge (SOC) as shown in Figure 4. However, if the battery state of charge is lower than a minimum threshold (SOCmin), the battery charge is ensured at the nominal charge current (Ich_nom). In the case where the battery state of charge exceeds a maximum threshold (SOCmax), the battery charge is disabled.

The block diagram of AFEM is given in Figure 5. A low pass filter is used to determine the low frequency part of the mission. Practically, in order to overcome discontinuity problems in the FC mission (I_FC), two cutoff frequencies are attribuated to the low-pass filter. The first frequency (Fg1) is used during the positive current gradient phases. The second (Fg2) is only used during the negative current gradient phases. It should be noted that the Fg2 frequency should be higher than the first frequency (Fg1) so that the output of the low pass filter nearly follows the initial mission (I_load) during the negative current gradient phases, as quickly as possible.

Thus, as illustrated on the synoptic of Figure 5 for AFEM strategy, a first part of the FC mission (I_FC0) is equal to the low frequency part of the load mission (I_LF) obtained with the frequency switching between (Fg1) and (Fg2) following the sign of the load current gradient. A second part of the FC mission is related to the ''Battery current charge loop'' block. This latter block allows maintaining the battery SOC despite the battery system losses and the asymmetry of the mission profile. This same battery current charge (Ich_BAT) is deducted from the initial battery mission (I_HF).

Note that in the case of symmetric energy management strategy, the unique cut-off frequency Fg represents a degree of freedom (both the FC profile and battery profile are sensitive to this frequency) which affects both fuel cell and storage sizing. In the context of an integrated optimal design, this frequency is considered as an optimization parameter.

In the case of asymmetric energy management strategy, the cut-off frequency Fg1 is determined in the same way (by optimization) as the cut-off frequency Fg in the case of symmetric energy management strategy. It is necessary that the cut-off frequency Fg2 is high simplify the analysis, the case considered is where the battery state of charge is higher than the maximum threshold (SOCmax). The battery current charge loop is then disabled

(Ich_BAT = 0).
It should be noted that only a simplified system vision has been put forward in this section in order to emphasize the comparative performance of both (SFEM vs AFEM) management strategies. Thus, current controllers have not been described as being composed of classical structures. In reality, these controllers are necessary to control currents and also to limit current peaks during transients. However, these constraints are not so drastic due to the relatively slow dynamics of the practical mission profile. Furthermore, these control loops have been implemented in experiments discussed in section 4.

Finally, as for the SFEM, it is also true that regenerative currents would set constraints on battery charge currents even with the AFEM, leading to a storage oversizing. However, as for our case study, except to the hybrid vehicle case, a wide class of embedded systems does not require managing regenerative powers; it is for example the case of aircraft systems for which the load profile does not contain energy recovery phases (negative charge power). This application field is in fact our main target but this issue has not been emphasized due to generality concern of our talk.

The simulation of the hybrid system

In this section, the simulation results of the hybrid system with both energy management strategies are presented. The simulation software used in this study is Matlab Simulink.

The hybrid system simulation is based on quasi-static electrical models as illustrated on Figure 6. According to our scientific experience [START_REF] Rafal | Experimental Validation of Hybrid Emergency Electrical Network for Aircraft Application[END_REF], [START_REF] Roboam | Systemic Design Methodologies for Electrical Energy Systems: Analysis, synthesis and management[END_REF], the global energy behavior at system level is quasi equivalent whatever dynamic electrical models with control loops or quasi-static electrical models with ideal current controllers are used. This remark is valid only if current controller bandwidths are rapid enough to deal with the reference dynamics, which is the case in our study due to the medium dynamics on the mission profile. As an example, dynamics displayed in Figure 3 (i.e some seconds of time range) to illustrate power management issues are quite slow with respect to current controller time responses (i.e some milliseconds of the time range). This model simplification is also validated through our experimental study (Figure 11) which suggests the same behavior of the hybrid system in terms of current sharing between battery and fuel cell devices. Experimentally, we have also obtained a good bus voltage quality and stability (Figure 12) which proves that voltage/current controllers are sufficiently fast to be idealized at the hybrid system level, making relevant the assumption of quasi-static electrical models.

The fuel cell quasi-static electrical model takes into account the activation (ηact), ohmic (Relec×IFC) and diffusion (ηdiff) losses. The fuel cell voltage is obtained by subtracting the total losses from the reversible potential Erev as shown in Figure 6.a.

The input of this model is the current I_FC stemming from energy management strategies given in Figure 2 or Figure 5.

Concerning the battery quasi-static electrical model, a R, E electrical model is used to obtain the battery cell voltage vBAT. The input of this model is the current I_BAT stemming from energy management strategies given in Figure 2 or Figure 5. The cell current is obtained by dividing I_BAT by the number of battery blocks in parallel (NPBAT).

The battery resistance rBAT and its electromotive force eBAT are interpolated from manufacturer's data according to the cell state of charge (SOC). The battery SOC is obtained from the ratio of the cell charge quantity (qBAT) to its nominal capacity (C). The charge quantity (qBAT) is defined as the sum of the integral of the cell current (iBAT) and the initial charge quantity of the cell (q0). According to the sign convention, the battery discharge current is considered to be positive.

Ohmic losses

I_FC I_FC FC elec diff act rev FC I R E V       + -

Diffusion losses Activation losses

Reversible potential The simulated profile (I_load) corresponds to the actual load profile to be fulfilled by the hybrid system. A cycle of 12 s (Figure 7) is continously repeated.

E rev η act η diff R elec v BAT I_BAT 1 BT v BAT r BAT BAT i r BAT e  BAT q C 1 soc 0 q -1 ) (soc e BAT ) (soc r BAT N PBAT i BAT
Figure 7.a shows simulation results of the hybrid system with SFEM. The cutoff frequency is fixed at 0.3 Hz. This value is compatible with the response-time of an H2/O2 FC system selected for the tests in question. Verification was carried out to ensure that the FC provides the low frequency part of the mission while the battery provides the rest (i.e the high frequency part).

It should be remembered that in the present study, the hybrid system topology shown in Figure 8 was initially sized. This system is composed of a 5 kW H2/O2 FC stack with a Li-Ion battery characterized by a nominal voltage of 216 V and a nominal capacity of 6 Ah (the series connection of 60 accumulator cells with a nominal voltage of 3.6 V). As recommended by the manufacturer (due to ageing issues), the battery may be solicited during discharge with a maximum current of 15 C (90 A) and with a maximum charge current of 1 C (6 A). Under such conditions, the 216 V battery is unable to satisfy the current profile I_BAT shown in Figure 7.a. Indeed, the charge current (the negative part of the curve) strongly exceeds, in absolute value, the nominal charge current (6 A). In order to satisfy this frequency sharing of the mission, one solution is to use two 216 V batteries in parallel. This solution obviously leads to a significant weight increase.

Unlike SFEM, the mission sharing according to AFEM fully respects the source sizing. In fact, Figure 7.b shows that the battery charge current (the negative part of I_BAT) does not exceed the nominal charge current (Ich_nom = 6 A). It can be observed that the level of the battery charge current depends only on the battery state of charge. Some of the dynamic constraints (strong negative current gradients for the I_FC in Figure 7.b) are transferred on the FC a priori without any impact on its lifetime. Further aging tests are forecast to confirm this assessment.

We verify in Figure 7.b that the FC guarantees the low frequency part of the mission only during the positive current gradient phases. It also ensures the battery charge current.

Experimental results

As previously mentioned, the LAPLACE laboratory is equipped with a ''Hydrogen'' platform with facilities used here for testing a 5 kW H2/O2 FC stack with a 216 V Li-Ion battery. Both sources are connected to a 270 V HVDC bus through two DC/DC power converters (Figure 8); it should be noted that the 270 V DC bus is called "HVDC" in the aeronautic field (which is our main application target) contrary to low voltage DC for 28 V avionic loads. The test bench structure ''FC-Li-ion Battery Hybridization'' is illustrated in Figure 8 with some pictures of devices in Figure 9. Figure 10 shows also the evolution of the battery SOC. In the first 5 minutes, the battery SOC decreases rapidly from 100% to 81%. In fact, during this period, the battery is the only available power source of the hybrid system. The battery recharging loop is not yet active. At t = 5 minutes, the hybridization starts and the battery recharging loop is activated which allows maintaining the battery SOC around 60%. At t = 2 hours and 42 minutes, the The losses in the fuel cell power converter are compensated by the fuel cell by producing an additional power resulting in H2 and O2 over-consumption. Concerning the battery system losses (battery losses + battery converter losses), a loss compensation loop is required to prevent the drift of the battery SOC. These losses are also supplied by the fuel cell. In the case of SFEM, the dynamic of the loss compensation loop should be slow regarding the fuel cell dynamic. In the case of AFEM, the loss compensation loop is integrated in the battery recharging loop. In both cases, the compensation of the battery system losses leads to H2 and O2 over-consumption.

The waveforms of both "HVDC bus side" FC and battery currents, respectively referred to as I_FC_HVDC and I_BAT_HVDC, and the load current I_Load at the HVDC bus are given in Figure 10 and Figure 11. The experimental results are perfectly in accordance with AFEM simulation results given also in Figure 10 and Figure 12 shows the experimental evolution of the DC bus voltage. This result shows a good bus voltage stability. This latter issue shows that voltage/current controller bandwidths are fast enough to face the load profile dynamics. In that particular case, the bus voltage has been regulated by the battery storage sub system while the FC sub system is current controlled by fulfilling the current sharing strategy as proposed in AFEM principle. In this section, the optimal design of the hybrid system is discussed; the optimized process is the one displayed in Figure 8. The objective of this design approach is to identify solutions with low weight and low losses.

The optimization process

The optimization process is given in Figure 13. A Non-dominated Sorting Genetic Algorithm (NSGA-II) [START_REF] Deb | A fast-elitist non-dominated sorting genetic algorithm for multiobjective optimization: NSGA-II[END_REF], [START_REF] Sareni | Fitness sharing and niching methods revisited[END_REF] is used to solve this two objective problem. This algorithm allows generating the design variables which are then injected into the hybrid system sizing model in order to evaluate the optimization criteria and constraints. The load current profile is an input for the simulation of the hybrid system model.

The design variables

The design variables (here 9 variables) are illustrated in Figure 14 and are given in Table 1. With regard to the energy source sizing models (FC and battery), the reader is invited to consult [START_REF] Chanpeng | Model and experiment analysis of 1.2 kW PEMFC electrification[END_REF], [START_REF] Arregui | Theoretical study of a power generation unit based on the hybridization of a fuel cell stack and ultracapacitors[END_REF], [START_REF] Jaafar | Sizing and energy management of a hybrid locomotive based on flywheel and accumulators[END_REF] and [START_REF] Menard | Dynamic Energy Model of a Lithium-Ion battery[END_REF]. The sizing model of the power converter is given in [START_REF] Mariani | Design of a medium voltage power converter storage devices embedded in a hybrid emergency network for more electrical aircraft[END_REF].

The design constraints

The optimization constraints are related to technological limitations of energy sources and to the satisfaction of the current sizing profile. Five inequality constraints formulated in terms of minimization (i.e. gi ≤ 0, i = {1,…,5}) should be respected to guarantee the feasibility of the hybrid system.

As the number of battery cells is the product of two design variables (NBAT = NSBAT × NPBAT), the equivalent battery should provide a suitable response to some energy and power constraints:

 The maximum discharge power

The maximum discharge power that the battery can provide (Pdisch_BAT_max) should be higher than or equal to the maximum discharge power of the battery profile (Pdisch_profileBAT_max).

g1 = Pdisch_profileBAT_max -Pdisch_BAT_max ≤ 0 (1)
 The maximum charge power

The maximum charge power authorized by the battery (Pch_BAT_max) should, in absolute value, be higher than or equal to the maximum charge power of the battery profile (Pch_profileBAT_max).

g2 = │Pch_profileBAT_max│ -│Pch_BAT_max│ ≤ 0 (2)
 The stored energy

The energy which can be stored in the battery (Etot) should be higher than the energy required for the battery profile (Eu). İndeed, according to the manufacturer's recommendations, only a percentage of Etot is allowed to satisfy a certain maximum Depth of Discharge (DODmax) [START_REF] Belouda | Integrated optimal design and sensitivity analysis of a stand alone wind turbine system with storage for rural electrification[END_REF]:

g3 = Eu -DODmax×Etot ≤ 0 (3)
Concerning the FC stacks, the constraints are:

 The number of FC cells in series

For structural reasons and gas distribution problems, the number of FC cells in series is limited to 300 [36], [37]:

g4 = NFC -300 ≤ 0 (4)

 The FC surface

The FC surface is also subjected to technological constraints related to current density, uniform distribution of temperature and mechanical clamping forces. It is limited to 2500 cm 2 : g5 = SFC -2500 ≤ 0 (5)

The design objectives

Two conflicting design objectives are considered in this study: the system weight F1 and the energy of system losses F2.

   CV-BAT BAT CV-FC FC H CV-BAT BAT CV-FC FC + Losses + Losses + Losses = Losses F + M + M + M + M = M F 2 2 1 (6)
where MFC, MCV-FC, MBAT, MCV-BAT respectively denote FC stack mass, FC power converter mass, battery mass and battery power converter mass. MH2 denotes the mass of hydrogen consumed during the mission and the mass of the corresponding tank. The stored hydrogen density is fixed to 5%: a 100 kg tank is needed to store 5 kg of H2 [START_REF] Adams | Identification of the optimum on-board storage pressure for gaseous hydrogen city buses[END_REF].

The total losses energy F2 is obtained by the time integral of the power loss in both energy sources (FC stack and battery) and power converters.

The optimization results

A comparison between the optimal sizing of the hybrid system with SFEM and that which is obtained with AFEM has been carried out. Both Pareto fronts are given in Figure 15. For confidentiality reasons, the results of Figure 15 are given in per-unit. It is clear that asymmetric frequency energy management (AFEM) strongly dominates symetric strategy (SFEM); it offers significant gain in terms of system weight compared to symmetric frequency energy management (SFEM).

Figure 16 shows the distribution of the system mass between the battery system (battery converter + battery) and the FC system (FC stack + FC converter + FC auxiliaries + H2 tank). This difference between the two management strategies is mainly due to the batteries whose sizing constraints, in the case of SFEM, are due to the charge current: as recommended by the manufacturer the maximum battery charge current is equal to 1 C (6 A) to preserve life duration. Indeed, the battery mass represents approximately 70% of the system mass. The FC system mass is almost the same with the two energy management strategies (AFEM and SFEM).

It should be noted that in order to minimize the system losses, the optimizer tends to choose oversized fuel cell stacks. Indeed, from the energy point of view, a large-sized FC is more efficient than a small one.

By making the assumption that the battery can accept higher charge currents, a significant gain in terms of mass may be achieved compared to a charge current of only 1C. However, responding to fast charging phases may possibly deteriorate the battery lifetime. SFEM: Battery system mass SFEM: FC system mass AFEM: Battery system mass AFEM: FC system mass Fig. 16. The system mass distribution between the battery system (battery converter + battery) and the FC system (FC stack + FC converter + FC auxiliaries + H2 tank)

Conclusion

In the present paper, two energy management strategies of an autonomous hybrid system based on the association of a Fuel Cell and a Li-Ion battery are presented. Both strategies are based on a frequency sharing of the mission; here the current profile which is related to the power sharing given the bus voltage. This frequency based strategy ensures compatibility between frequency components of the mission and intrinsic dynamical characteristics of sources, especially when the power profile is positive without regenerative phases; this kind of profile is usually required in various classes of embedded systems as in aircraft applications. Through simulations of a real (typical) mission of the hybrid system, the authors have shown the interest of Asymmetric Frequency Energy Management compared to Symmetric Frequency Energy Management in terms of system sizing. A number of dynamical constraints is transferred to the FC; conclusive evidence related to the lack of lifetime impact

is not yet available. As a following, experimental results of Asymmetric Frequency Energy Management are presented. Finally, the impact of the choice related to the energy management strategy on the hybrid system weight has been illustrated in the context of an integrated optimal design: this study confirms the interest of an original asymmetric strategy for this example of typical source association and in the framework of typical mission profiles.

Figure 2

 2 Figure 2 illustrates the principle of SFEM. The variables I_BAT and I_FC represent the mission parts respectively sent to the battery and to the FC. Wg denotes the cutoff pulsation of the low pass filter (Wg = 2πFg).

Figure 3 .Fig. 2 .Fig. 1 .

 321 Figure 3.a shows an example of frequency sharing applied to an "academic" mission (square wave) with a cutoff frequency Fg = 0.5 Hz. The FC provides the low frequency part of

Fig. 3 .

 3 Fig. 3.b. Illustration of Asymmetric FrequencyEnergy Management

Fig. 4 .

 4 Fig. 4. The principle of the battery recharging loop

Figure 3 .

 3 Figure 3.b illustrates the principle of AFEM with a cutoff frequency Fg1 = 0.5 Hz. To

Fig. 5 .

 5 Fig. 5. Principle of the Asymmetric Frequency Energy Management

Fig. 6 .

 6 Fig. 6.a. Fuel cell quasi-static electrical model Fig. 6.b. Battery cell quasi-static electrical model

Fig. 8 .

 8 Fig. 7.b Simulation of the mission profile with Asymmetric Frequency Energy Management

Fig. 9 .

 9 Fig. 9. The ''FC-Li-ion Battery Hybridization'' test bench on the « Hydrogen » platform of the LAPLACE laboratory

Figure 11 Fig. 10 .Fig. 11 .

 111011 Fig. 10. The experimental and simulation results of Asymmetric Frequency Energy Management strategy

Fig. 12 . 5 .

 125 Fig. 12. The experimental evolution of the DC bus voltage

Fig. 13 .

 13 Fig. 13. Optimization process

Fig. 14 .

 14 Fig. 14. The illustration of design variables

Fig. 15 .

 15 Fig. 15. Comparison of the two Pareto fronts obtained with symmetric and asymmetric energy management strategies

Table 1 .

 1 List of design variables

	Symbol	Design Variable
	Fg	Power sharing frequency (the cutoff frequency of the low pass filter)
		Fg = Fg1 in the case of AFEM (Fg2 being fixed at 100 Hz)
	NSBAT	Number of battery cells in series
	NPBAT	Number of battery blocks in parallel
	Nb//-BAT	Number of interleaved battery converter branches
	Fd-BAT	Battery converter switching frequency
	Vstack-min	FC stack minimum voltage
	JFC-max	FC maximum current density
	Fd-FC	FC stack converter switching frequency
	Nb//-FC	Number of interleaved FC stack converter branches


F Sambou

J J Gonzalez

M Benmouffok

P Freton

Theoretical study of the arc motion in the hollow cathode of a dc thermal plasma torch

Keywords: hollow cathode, arc motion, 3D model

A 3D model of a hollow cathode of a thermal plasma torch was developed to study the motion of the electric arc in the upstream part as well as its influence on the characteristics of the thermal plasma generated at the entrance of the downstream part (anode). The translational motion along the axis and the rotational velocity of the cathodic arc root are studied. The results of the model are presented and compared with those of the literature and also with experimental results. Good agreements on the prediction of the axial arc position and of the rotational speed of the arc root are found with the literature. A correlation is also shown between the theoretical positions of the arc root and the experimental traces of erosion observed on the cathode.

Introduction

Thermal plasma torches with hollow cathode are used in several industrial applications such as: waste treatment, production of reducing gases in the steel industry, petro chemistry [START_REF] Boulos | Springer Science, Fundamentals and applications[END_REF][START_REF] Fauchais | [END_REF][START_REF] Laroche | Book of Les plasmas dans l'industrie[END_REF]. This plasma torch technology is used for high power outputs in the MW range [START_REF] Brilhac | [END_REF] and allows the generation of large volumes of thermal plasma at very high enthalpies in the range 107-108J/kg and temperatures between 3000K and 15000K [START_REF] Boulos | Springer Science, Fundamentals and applications[END_REF]. These torches have some specificities that make their study delicate: (1) a hollow cathode, (2) a vortex gas injection between the electrodes, (3) an external magnetic field in order to rotate the arc root in the upstream electrode and to limit its erosion. The complexity of their geometry and the extreme conditions prevailing in these torches (High Power) make their experimental studies difficult and expensive [START_REF] Brilhac | [END_REF]. The experimental characterizations of hollow torch technologies are devoted to measurements of macroscopic physical quantities such as: arc voltage, emitted sound waves, and temperature and velocity fields of the jet at the torch exit. It does not allow a detailed understanding of their internal functioning as the flow path from the injection chamber or description of the arc attachment in the hollow cathode. In the anode region, the difficulty to visualize the arc anode root or to quantify the erosion is also difficult, however studies performed not in hollow cathode but in rod cathode have being published showing new interesting results and possibilities [5,6]. Because of these limitations in hollow cathode, numerical studies are the most suitable for a better understanding of their internal behavior.

In the literature, few studies on this torch technology exist probably due to the complexity of their configuration. In the 1990s, Brilhac et al [START_REF] Brilhac | [END_REF] carried out a theoretical and experimental study of a thermal plasma torch with a hollow cathode for a power less than 100 kW and without an external magnetic field. During their experimental study, they noted the positioning of the wear zones as a function of the operating conditions and linked the position of this zone to the variations of the axial component of the velocity in close proximity to the cathode surface (0.5mm from the surface) using a cold theoretical analysis of the plasma gas flow (in the absence of electric arc) using the numerical code "MELODIE". Thanks to their experimental electrical and acoustic measurements, they qualitatively stipulate that the variations in the output characteristics of the torch are related to the fluctuations of the arc within the cathode. In the 2000s, again through the experiment, Hur et al [7][8] studied the behavior of a 200kW thermal plasma torch in the absence of an external magnetic field. Based on the same principle as Brilhac et al [START_REF] Brilhac | [END_REF] concerning the analysis of the cold plasma gas flow and thanks to their experimental measurements they also correlated the position of the wear zone to the variations of the axial component of the cold plasma gas flow velocity near the cathode. This qualitative technique for determining the position of the wear zone on the surface of the cathode by Brilhac et al [START_REF] Brilhac | [END_REF] is also used in configurations in the presence of a magnetic field by Park et al [9]. In their study, they consider the external magnetic field to be axial and constant. In the 1990s Minoo et al [10] conducted an experimental and theoretical study on this type of torch technology and showed that the positioning of the arc root on the cathode surface (position of the wear zone) results from an equilibrium between the magnetic and drag forces of the arc root by the flowing plasma gas. Another study in 2009 by Freton et al [11] shows that the external magnetic field, in particular its radial component, has a strong influence on the axial position of the arc root on the cathode surface. They find through their results that the external magnetic field is far from constant and that the arc root tends to reach a position where the radial component of the external magnetic field is at its maximum.

In this paper, we present a 3D numerical model that we have developed to represent the motion of the arc and its influence on the thermal plasma characteristics. We focus the study only on the upstream part of the torch (cathode). The first part of this paper will be devoted to the presentation of the numerical model. In the second part, a particular attention will be devoted to the calculation of the magnetic field generated by the external coil. The third part will present and discuss the results of the model for given operating conditions with a comparison with the literature. Finally, in the last part, a conclusion will be given.

Numerical model

In the framework of this study, the numerical tool used to solve the equations is the commercial code Ansys Fluent, which is widely used in the literature for the study of fluid flows at low temperatures, by solving the Navier-Stokes equations of fluid mechanics. This software has functionalities allowing the user to add additional terms to the Navier-Stokes equations (source terms) as well as additional equations. These functionalities, called UDF (User Define Function), are external functions that the user will be able to define and implement in the software to adapt it to his problematic. Thus in the framework of our study, these functions are used for the description of plasmas properties, specific source terms and for the electromagnetic aspect of thermal plasmas. To solve the equations, this tool uses the finite volume method proposed by Patankar et al [START_REF] Patankar | Book of Numeritical Heat transfer and Fluid Flow[END_REF]. According to this method, the equations are rewritten in the form of the generalized differential equation (1): 𝜕(𝜌𝜙) 𝜕𝑡 + 𝑑𝑖𝑣(𝜌𝑣 ⃗⃗⃗ 𝜙) = 𝑑𝑖𝑣 (Г 𝜙 . 𝑔𝑟𝑎𝑑 ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ (𝜙)) + 𝑆 𝜙 [START_REF] Boulos | Springer Science, Fundamentals and applications[END_REF] In this equation, we have respectively the temporal term, the convective term, the diffusive term and the source term ( 𝑆 𝜙 ). 𝜙 is the physical property studied, Г 𝜙 its diffusion coefficient, 𝜌 the mass density, and 𝑣 ⃗⃗⃗ the velocity.

Hypothesis and equations

2.1.1. Hypothesis. We use air as plasma gas. We consider that the thermal plasma created is a Newtonian fluid in laminar flow, at local thermodynamic equilibrium (LTE). The LTE hypothesis is used for the calculation of the properties of the thermal air plasma [START_REF] Bacri | [END_REF] such as thermodynamics, transport and radiative properties. This assumption is classically adopted in the plasma core, nevertheless departures from equilibrium can exist in plasma surrounding and in the electrodes sheathes. According to P. Andre et al [X], for heavy species temperature T=10kK, the frequency collision of the electron with other particles is at equilibrium equals to 10 13 s -1 and equals to 5.10 13 s -1 for =Te/Th=3 (ratio of electronic temperature on heavy one). This number of collisions by second is important compared to the calculation time classically used in the simulation (10s in our case) allowing considering LTE assumption. For validating the laminar flow hypothesis, we use the method suggested by Sinkevich et al [START_REF] Sinkevich | Boundaries of stability of the plasma flow in cylindrical channels, TPP 7 Thermal Plasma Processes[END_REF] for thermal plasma to estimate the Reynolds number. Indeed, in the case of thermal plasma, the transport and thermal proprieties change highlight with the plasma flow dimensional characteristics (essentially in the radius direction for a tubular channel). This highlight change makes difficulty to choose the characteristics of the plasma flow allowing the Reynolds number estimation. Therefore, this method suggests estimating the Reynolds number by calculating mean values of the characteristics of the plasma flow. According to this method, we obtain a Reynolds numbers within the torch below 500. That allows us to consider a laminar flow hypothesis in our study. The radiation source term is taken into account as a loss term in the energy equation by the net emission coefficient method [START_REF] Naghizadeh-Kashani | [END_REF]. The walls erosion due to the arc attachment or incoming radiation from the arc is not taken into account.

Equations. The equations solved with the source terms are as follows:

© xxxx IOP Publishing Ltd -The mass conservation equation :

∂ρ ∂t + ∇ ⃗ ⃗ . (ρ𝑣 ) = 0 (2)
Where t represents the time, ∇ ⃗ ⃗ is the spatial gradient operator, 𝜌 is the mass density, and 𝑣 ⃗⃗⃗ is the vector velocity -The momentum conservation equation:

∂(𝜌𝑣 ⃗ ) ∂t + ∇ ⃗ ⃗ • (𝜌𝑣 𝑣 ) = ∇ ⃗ ⃗ (𝜇∇ ⃗ ⃗ • 𝑣 ) + 𝑆 𝑣 (3)
Where

𝑆 𝑣 = -∇ ⃗ ⃗ (𝑃) + ∇ ⃗ ⃗ • 𝜏̿ + 𝑗 × 𝐵 ⃗ + 𝜌𝑔 ( 4
)
𝑆 𝑣 is the momentum source term where we distinguish respectively pressure forces (𝑃), viscous stresses 𝜏 ̿, Lorentz forces (product of the current density by the magnetic field 𝑗 × 𝐵 ⃗ ) and gravity 𝑔 . 𝜇 being the viscosity. The stress tensor 𝜏̿ is given for a Newtonian fluid by the relation:

𝜏 𝑖𝑗 = 𝜇 [( ∂𝑣 𝑖 ∂j + ∂𝑣 𝑗 ∂i ) - 2 3 ∂𝑣 𝑖 ∂i 𝛿 𝑖𝑗 ] (5)
Where the index 𝑖 and 𝑗 correspond to the different components of the coordinate system, 𝑣 𝑖 and 𝑣 𝑗 are the velocity components. The symbol 𝛿 𝑖𝑗 = 1 if 𝑖=𝑗 and 𝛿 𝑖𝑗 = 0 elsewhere -The energy conservation equation:

∂(ρh) ∂t + ∇ ⃗⃗⃗ • (ρ𝑣 ⃗⃗⃗ ℎ ) = ∇ ⃗⃗⃗ ( К 𝐶 𝑝 ∇ ⃗⃗⃗ ℎ) + 𝑆 𝐻 (6)
Where

𝑆 𝐻 = 𝑗 ⃗⃗ 2 𝜎 -𝑈 + 5 2 𝑘 𝐵 𝑒 𝑗 ⃗ ⃗ • ∇ ⃗⃗⃗ ( 𝑇 𝐶 𝑝 ) ( 7
)
𝑆 𝐻 is the energy source term, where we have respectively the joule effect (Ratio of the square of current density on the electrical conductivity 𝜎), radiative losses (𝑈), and enthalpic flux of the electrons. 𝑘 𝐵 is the Boltzmann constant, 𝐶 𝑝 the specific heat at constant pressure, Ƙ the thermal conductivity, ℎ the enthalpy, 𝑇 the temperature and 𝑒 the elementary charge.

According to the net emission coefficient method, the radiative losses are given by:

𝑈 = 4𝜋𝜀 𝑁 (8)
Where 𝜀 𝑁 represents the net emission coefficient. For the description of the electromagnetic phenomena of the plasma, it is necessary, in addition to the fluid equations, to solve the equation of the electric potential and the vector potential equations. These equations allow the calculation of electric and magnetic fields respectively. These electromagnetic fields act on the charge carriers of the plasma (electrons and ions) and consequently lead to the coupling between the fluid and electromagnetic equations by acting as sources (through the Lorentz force) in the momentum equation (3) and energy (joule effect) in the energy equation (6). The electromagnetic equations are established by assuming that the plasma medium is a succession of stationary states.

In order to satisfy the current continuity equation, we use the equation (9) where V designates the scalar potential.

𝛻 ⃗ • (𝜎𝛻 ⃗ 𝑉) = 0 (9)
Equation ( 9) is used to calculate the electric field (10) and the current density (11):

𝐸 ⃗ = -𝛻 ⃗ 𝑉 (10)
𝑗 = 𝜎𝐸 ⃗ (11)
© xxxx IOP Publishing Ltd

For the calculation of the magnetic field 𝐵 ⃗⃗⃗ , two methods are generally used in the literature [16]:

-Solving the potential vector equation 𝐴 ⃗⃗⃗ using the Laplacian ∆𝐴 :

∆𝐴 = -𝜇 0 𝑗 ( 12
)
𝜇 0 is the magnetic permeability of vacuum

-Biot & Savart law : 𝐴 (𝑟 ) = 𝜇 0 4𝜋 ∭ 𝑗 (𝑟 ′) |𝑟 -𝑟 ′| 𝑑𝑟′ 3 ( 13
)
Where 𝑟 is the position of the point where we want to calculate the magnetic field and 𝑟 ′ is the position of the current element contributing to this magnetic field. The magnetic field is deduced by the vector potential according to the relation ( 14):

𝐵 ⃗ = 𝛻 ⃗ × 𝐴 (14)
In this study, we use a Cartesian coordinate system. The set of conservations equations solved by the model are the same equations described by Freton et al [11].

Geometry

In the figure 1, we present the cathode geometry of the studied torch. The torch is made up of two electrodes (anode and cathode) and an injection chamber composed by 6 injectors inclined at an angle of 20° through which the plasma gas is injected. The anode and the cathode are electrically isolated by insulators in the injection chamber. This configuration of the injectors makes it possible to produce a vortex-type flow of the plasma gas. The torch is also equipped with an external coil in the cathode side. The internal details of the torch constrain the plasm formation, nevertheless the external magnetic field produced by the coil and the localization of its maximum value depending on the coil position act on the arc root rotation and stabilization. The management of the downstream arc over the anode can also in some configurations be important

© xxxx IOP Publishing Ltd
but in our studied configuration it doesn't exist in the anode side external magnetic field produced by a coil. The cathode's coil is represented in the Figure 1. Note that for the geometry studied, the anode part has been truncated just after the injection chamber in order to focus the study of the arc motion only in the cathode. Gris parts represent the coil, and the torch bottom is located in the right side. In this study, the torch dimensions are expressed in adimensional length unit (a.l.u) as preconized by our industrial partner and for confidentialities reasons. However the geometry is drawn as per scale maintaining the aspect ratio. The total length of the geometry presented in Figure 1 is about 30cm. The inlet diameters of the hollow cathode and of the anode are respectively dc =0.1 a.l.u., da= 0.083 a.l.u and dch = 0.26 a.l.u (Figure 1). 1 presents the boundary conditions.

Boundary conditions. Table

Table 1. Boundary conditions

Boundary Injectors inlet Torch exit Walls T[K]

300

𝛻 ⃗ 𝑇 • 𝑛 ⃗ = 0 300 𝑫 [g.s -1 ] 37 - - 𝑷 [atm] - 𝛻 ⃗ 𝑃 • 𝑛 ⃗ = 0 𝛻 ⃗ 𝑃 • 𝑛 ⃗ = 0 𝑽 [𝑽] 𝛻 ⃗ 𝑉 • 𝑛 ⃗ = 0 0 𝛻 ⃗ 𝑉 • 𝑛 ⃗ = 0 𝒗 𝒊=𝒙,𝒚,𝒛 [m/s] 𝛻 ⃗ 𝑣 𝑖 • 𝑛 ⃗ = 0 𝛻 ⃗ 𝑣 𝑖 • 𝑛 ⃗ = 0 𝛻 ⃗ 𝑣 𝑖 • 𝑛 ⃗ = 0 𝑨 𝒊=𝒙,𝒚,𝒛 [T.m] 𝛻 ⃗ 𝐴 𝑖 • 𝑛 ⃗ = 0 𝛻 ⃗ 𝐴 𝑖 • 𝑛 ⃗ = 0 𝛻 ⃗ 𝐴 𝑖 • 𝑛 ⃗ = 0
In Table 1, D is the inlet mass flow rate and 𝑣 𝑖=𝑥,𝑦,𝑧 and 𝐴 𝑖=𝑥,𝑦,𝑧 the three components respectively of the velocity and vector potential. The current, after passing in the external coil, flows through the torch in the solid part. Depending on the plasma conditions in the fluid part, the current path is self-determined. We therefore apply an arbitrary current density boundary condition in the solid part of the cathode (𝐽 crown) on a crown located below the bottom of the cathode (figure 2).

On the downstream section (Torch exit, Figure 1) of the anode portion, we consider a porous anode condition [11] that is to say a zero electric potential boundary condition where the arc column will hang (Figure 2). In order to take into account the contribution of the magnetic field due to the current flow in the torch frame, Biot and Savart's law (equation 13) was used as a closing boundary condition on the solid boundaries for the calculation of the vector potential. In the solid, we solve only the scalar potential. In the Figure 2, we present a general view of the mesh. The mesh has been refined at the electrodes surfaces (boundary between fluid and solid parts) as well as at the injectors region in order to correctly capture the physical phenomena in these areas. To take into account the current flow in the cathode, the solid part of the cathode is also meshed. After a specific study on the mesh to obtain an invariant solution on converged results, the total resulting mesh has about 670 000 cells with a volume range from 6.3.10 -2 mm 3 to 3.4.10 -1 mm 3 . The larger cells are located at the top of the hollow cathode far from the arc presence as the smaller cells are in the plasma core on the walls where the arc root is attached. It can be seen in Figure 2 that the mesh is refined in the plasma area where will be located the plasma. It consists of two types of mesh: (1) an un-structured mesh: in the injection chamber and at the cathode extremity; (2) a structured mesh: in the electrodes as well as in the inter-electrode space.

2.2.2.

Procedure for initializing the calculation. In order to allow the flow of electric current in the plasma gas between the two electrodes for arcing, a hot conductive channel must be defined. This is done by introducing one energy source term via a UDF. The conductive channel was ignited between the axial positions z= 0.67 a.u.l. and z= 1 a.u.l. from the bottom of the cathode localized at z= 0 a.u.l., taken as the origin of the space in the model. The volume of conducting channel is approximatively 3.17 10 -6 m 3 and the instantaneous power deposited in the initial channel is around 125kW. This channel corresponds to an arc length and to one diameter equal respectively to 0.1m and 6.10 -03 m. The power is introduced as source term in the energy equation in all the cells volumes of the channel leading in our case to a mean initial temperature around 10kK. The initial temperature needs to be enough to allow a conducting channel. Regarding to the electrical conductivity variation versus temperature values from 9kK can be chosen. Higher values can be considered but the deposited energy must be negligible compared to the injected power in order to not introduce energy which could be longer to dissipate. Then the fluid and electromagnetic equations described above for the thermal plasma part are solved for 5 time steps (5 x 10𝜇𝑠) corresponding to a total deposited energy around 6.25J. Once the conductor channel has been established, we then run the whole calculation. The conductive channel obtained after the arcing is the one shown in Figure 3 and materialized by the 7kK temperature iso-surface. The temperature is in the range 7kK to 11kK. Of course the value of the deposited energy doesn't affect the results as after the channel establishment the power in the plasma is calculated by the imposed value of the current intensity and the resulting drop voltage directly calculated by the scalar potential resolution.

Arc attachment movements

In order to determine the position where the arc root is attached to the cathode surface, we consider the model proposed by Lindmayer et al [17] developed in the framework of the study of arc displacement in low voltage circuit breakers [18][19][START_REF] Lindmayer | IEEE Holm Conference on Electrical Contacts[END_REF][START_REF] Lindmayer | IEEE Holm Conference on Electrical Contacts[END_REF][START_REF] Ruimpler | 23rd International Conference on Electrical Contacts[END_REF][START_REF] Lindmayer | [END_REF][START_REF] Ruither | 23 rd International Conference on Electrical Contacts[END_REF]. This model is based on a calculation of the effective electrical conductivity at each point on the surface [17,[START_REF] Ms | PhD Thesis Modelization and Analysis of the Electric Arc in Low Voltage Circuit Breakers[END_REF][START_REF] Mutzke | 24 th International Conference on Electrical Contacts[END_REF] and allows introducing an additional sheath voltage drop the point of attachment having the highest conductivity. The rigorous numerical treatment of these zones would require the implementation of non-thermal equilibrium models that are more expensive in computing time [START_REF] Benilov | [END_REF] [START_REF] Cayla | Modélisation de l'interaction entre un arc électrique et une cathode[END_REF].

Once this threshold voltage has been reached and the arc has established itself between the electrodes, the voltage drop in these zones becomes independent of the increase in current density [17] and tends towards a constant voltage depending on the electrode material (10 Volt for copper). Thus, a thin layer of thickness ∆r is defined near the surface of the cathode. On this layer, a local value of the effective electrical conductivity 𝜎 𝑒𝑓𝑓 is calculated according to one of the voltage-current density characteristics [START_REF] Ms | PhD Thesis Modelization and Analysis of the Electric Arc in Low Voltage Circuit Breakers[END_REF] and the relation [START_REF] Naghizadeh-Kashani | [END_REF]:

𝜎 𝑒𝑓𝑓 (𝐽) = 𝐽 • ∆𝑟 𝑉 𝑠 (𝐽) (15)
The study of the motion of the arc by this model requires modeling the current path in the electrode, which imposes to mesh the electrode framework (solid part). The Figure 4 shows a schematic representation of the arc with the positions of the sheath. In [START_REF] Ms | PhD Thesis Modelization and Analysis of the Electric Arc in Low Voltage Circuit Breakers[END_REF] four curves (a, b, c and d) whose variations are dependent on the arc initiation threshold voltage (23.3V, 19.7V, 17.1V and 10V) can be used. The highest is the voltage peak; the most difficult is the motion of the arc. For the instantaneous displacement of the cathodic arc root, one considers at each instant through the model, the position of the arc root as corresponding to the weighted barycenter of each point of the cathode surface by its electrical conductivity, according to relations ( 16) and ( 17) respectively for the instantaneous axial and azimuthal positions :

© xxxx IOP Publishing Ltd

𝑍 𝑐𝑎𝑡 (𝑡) = ∑ 𝜎 𝑖 (𝑡)•𝑍 𝑖 (𝑡) 𝑛 𝑝 𝑖=1 ∑ 𝜎 𝑖 (𝑡) 𝑛 𝑝 𝑖=1 (16)
𝜃(𝑡) = ∑ 𝜎 𝑖 (𝑡)•𝜃 𝑖 (𝑡) 𝑛 𝑝 𝑖=1 ∑ 𝜎 𝑖 (𝑡) 𝑛 𝑝 𝑖=1 (17)
Where :

-𝜎 𝑖 (t) represents the electrical conductivity of point 𝑖 of the cathode surface;

-𝑍 𝑖 (𝑡) and 𝜃 𝑖 (𝑡) are respectively the axial and azimuthal positions of point 𝑖 on the cathode surface ; -𝑛 𝑝 is the total number of points on the cathode surface.

Calculation of the magnetic field generated by the coil

The external coil of this torch is a 3-layer solenoid type, each of the 3 layers having a different number of turns. To determine the distribution of the magnetic field generated by the coil we used an analytical method from the literature proposed by Conway [START_REF] Conway | [END_REF] adapted to our coil configuration. This method determines the magnetic field generated by a simple coil of finite dimension using elliptic integrals by integrating the Maxwell's equations of magneto statics, (18) and (19).

∇ ⃗ ⃗ × 𝐵 ⃗ = 𝜇 0 𝐽 (18)
∇ ⃗ ⃗ • 𝐵 ⃗ = 0 ( 19
)
𝐵 ⃗ represents the magnetic field, 𝐽 the current density, and 𝜇 0 the vacuum permittivity. By introducing relation [START_REF] Sinkevich | Boundaries of stability of the plasma flow in cylindrical channels, TPP 7 Thermal Plasma Processes[END_REF] in equation (19), we obtain the equation:

∇ ⃗ ⃗ • 𝐴 = 0 (20)
By introducing relation [START_REF] Sinkevich | Boundaries of stability of the plasma flow in cylindrical channels, TPP 7 Thermal Plasma Processes[END_REF] into equation ( 18) and using equation [START_REF] Lindmayer | IEEE Holm Conference on Electrical Contacts[END_REF], the vector potential satisfies to the potential vector equation [START_REF] Patankar | Book of Numeritical Heat transfer and Fluid Flow[END_REF]. In the case of an axisymmetric solenoid in cylindrical coordinates, only the azimuthal components of the vector potential (𝐴 𝜃 ) and current density (𝐽 𝜃 ) contribute to the generated magnetic field. Thus, equation ( 12) can be rewritten in the azimuthal direction by equation ( 21):

© xxxx IOP Publishing Ltd ∇ ⃗ ⃗ 2 (𝐴 𝜃 (𝑟, 𝑧)) = -𝜇 0 𝐽 𝜃 (𝑟, 𝑧) (21)
The relation ( 14) allows expressing the axial and radial components of the magnetic field as a function of the azimuthal component of the vector potential in cylindrical coordinates system ( 22) and ( 23):

𝐵 𝑧 (𝑟, 𝑧) = 1 𝑟 𝜕 𝜕𝑟 (𝑟𝐴 𝜃 (𝑟, 𝑧)) (22)
𝐵 𝑟 (𝑟, 𝑧) = -𝜕𝐴 𝜃 (𝑟,𝑧) 𝜕𝑧 [START_REF] Lindmayer | [END_REF] Note that there can be no toroidal component i.e. 𝐵 𝜃 (𝑟, 𝑧) = 0. By integrating equation ( 21) over the entire current distribution using elliptic integrals, Conway obtains the following analytical relationships for the axial and radial components of the magnetic field for a simple solenoid of finite dimension. These analytical relationships can be written by introducing the number of turns N of the solenoid by [11]:

𝐵 𝑟 (𝑟, 𝑧) = 𝜇 0 𝑁 𝐼 2 (𝐼 (0,1,1) (𝑅, 𝑟, |𝑍 2 -𝑧|) -𝐼 (0,1,1) (𝑅, 𝑟, |𝑍 1 -𝑧|)) (24)
𝐵 𝑧 (𝑟, 𝑧) = 𝜇 0 𝑁 𝐼 2 (𝐼 (0,1,0) (𝑅, 𝑟, |𝑍 1 -𝑧|) -𝐼 (0,1,0) (𝑅, 𝑟, |𝑍 2 -𝑧|)) if (𝑧 < 𝑍 1 ) (25)
𝐵 𝑧 (𝑟, 𝑧) = 𝜇 0 𝑁 𝐼 2 (2𝐼 (0,1,0) (𝑅, 𝑟, 0) -𝐼 (0,1,0) (𝑅, 𝑟, |𝑍 1 -𝑧|) -𝐼 (0,1,0) (𝑅, 𝑟, |𝑍 2 -𝑧|)) if (𝑍 1 < 𝑧 < 𝑍 2 ) (26) 𝐵 𝑧 (𝑟, 𝑧) = 𝜇 0 𝑁 𝐼 2 (𝐼 (0,1,0) (𝑅, 𝑟, |𝑍 2 -𝑧|) -𝐼 (0,1,0) (𝑅, 𝑟, |𝑍 1 -𝑧|)) if (𝑧 > 𝑍 2 ) (27)
𝐼 is the intensity of the current in the solenoid, R its radius, Z1 and Z2 the positions of its extremities. The functions 𝐼 (0,1,1) and 𝐼 (0,1,0) are the elliptic integrals or Bessel-Laplace's integrals [START_REF] Conway | [END_REF][START_REF] Abramovitz | Handbook of mathematical functions[END_REF].

Results

For this study, we applied an arc current intensity of 260A for an air flow rate imposed on the injectors of 37g/sec. The curve (c) [START_REF] Ms | PhD Thesis Modelization and Analysis of the Electric Arc in Low Voltage Circuit Breakers[END_REF] is used to consider the arc motion on a copper surface.

Magnetic field

In this section, we present the distribution of the components of the external magnetic field generated by the coil, calculated using the analytical relationships proposed by Conway et al [START_REF] Conway | [END_REF] adapted to our coil configuration. For this study, the axial component of the magnetic field responsible for the rotational motion of the arc is applied in such a way that it causes the cathodic arc root to rotate in the opposite direction of the vortex injecting plasma gas. Since the motion of the arc root is induced by the application of an external magnetic field, it is obvious that a reversal of the direction of the magnetic field should be able to impart a motion of the arc root just in the reverse direction [5]. Note that the total magnetic field within the torch represents the contributions of the magnetic field generated by the coil as well as the magnetic fields generated by the arc current and the cathode current.

Arc behaviour in the hollow cathode

In the Figure 7, we present the temporal behavior of the electric arc within the geometry studied, which is represented by the iso-surface 7kK. The sinuous plasma behavior part due to the external magnetic field, vortex injection, and changes of the arc root position is related to the hollow cathode (above red line) as the cylindrical part (bellow red line) is related to the portion of the anode part. This behavior of the electric arc is the consequence of the interactions between the blowing of the vortex airflow and its rotation under the action of the magnetic field. force diminishes and an axial equilibrium position exists after 10ms with a mean voltage value around Vmean = 626V, (Step 2) then the arc slowly oscillates from this axial position and still rotates due to the magnetic field. The rotational velocity keeps roughly the same, nevertheless the arc length has increased and the arc column submitted to convective and magnetic forces is more constraint leading to distortions and length variations illustrated by the amplitude of the voltage variations.

Cathodic arc root motion

In the Figure 9, we have represented on the left side the evolution of the axial position of the arc root on the cathode surface as a function of time and on the right side the axial profiles of the radial and axial components of the magnetic field generated by the coil. The position of the arc root at time t=0ms corresponds to its initialization position (z=0.67 a.l.u). The radial component of the external magnetic field being order smaller than the axial component, we had multiplied it by 10 in order to clearly see its spatial profile. This figure shows the instantaneous arc root axial position estimated by the equation 16. The arc root moves in a non-linear way towards the bottom of the cathode up to a mean equilibrium area where it oscillates. In this figure, the arc rotation equilibrium zone is delimited by the axial positons z=0.47 and 0.53 a.l.u. This non-linear axial movement of the arc root is the consequence of the arc fluctuations due to the joint actions of the plasma gas flow and the external magnetic field. The Figure 9 show also that the axial root attachment converges to a position located close to the maximum of the radial component of the magnetic field. Freton et al [11] have also observed a correlation between the rotational equilibrium area of the arc root with the position of the maximum values of the radial component of the external magnetic field.

In order to study whether there was a correlation between the erosion zones and the position of the rotational equilibrium zone, our industrial partner Europlasma carried out thickness measurements using an ultrasonic thickness gauge (from Krautkramer Branson). The measuring principle is as follows: Ultrasonic pulses are transmitted to the electrode material by a probe in contact with the surface of the cathode material. Some of the pulses are reflected at the surface in contact with the material while others pass through the material to its opposite side where they are reflected. By analyzing the paths taken by the sound pulses reflected at the contact surface of the probe with the material and those reflected at the opposite surface of the contact area, the thickness of the area of the material through which the ultrasonic pulses pass can be deduced. The electrode erosion is a slow process and the thickness of the eroded area changes versus time. This eroded area is so most visible after several hours of running. Nevertheless the eroded area position keeps the same. The accuracy in thickness measurements given in the technical documentation is about 0.25mm. In the figure 10 Figure 10 shows an overlap of the two zones: the theoretical hanging zone (arc root attachment: black frame 0.17 a.l.u. wide) and the experimentally observed one (wear zone: red frame 0.2 a.l.u. wide). Due to convective effects, the mean axial root position (Figure 9) is not located in the middle of the area of plasma interaction with the wall (Figure 10) as the plasma is blown by the upstream flow.

The arc interaction on the electrode leads to a heat transfer. Several physical quantities act as the radiation, the conduction and the electric components due to the current path between plasma and material. Depending on the arc motion and to the cooling system the wall temperature is more or less dissipated leading to erosion. In our model the sheath is not described and an "artificial method" is used to represent the sheath resistivity. Depending on voltage gradient [START_REF] Ms | PhD Thesis Modelization and Analysis of the Electric Arc in Low Voltage Circuit Breakers[END_REF] and on the curve chosen, more or less current is allowed to circulate through a calculated value of the electrical conductivity. The choice of the curve (a-b-c-d) could also allowed to adjust if necessary the dimension of the eroded zone. In case of departure from experimental results, the value of the ignition voltage has to be chosen on the basis of the measurement results [START_REF] Mutzke | IEEE Holm Conf. on Electrical Contacts[END_REF]. The higher is the ignition voltage, the smaller is the eroded area and the later begins the new arc root position. Nevertheless one can see in the figure 10 that our theoretical results on the prediction of the eroded zone show a similarity with the eroded zone found experimentally applying curve c. The figure shows a quasi-periodic azimuthal revolution of the arc around the axis of the torch. If the attachment is constricted there is a formation of arc root on the nozzle and the external magnetic field imparts a rotational movement to the constricted attachment due to appreciable Lorentz forces. This arc root movement is sometimes designed, induced root movement (IRM) [5]. Due to the forces the arc root rotates, each turn corresponds to one maximum in the Figure 11. The arc root pulls the plasma column which has a tortuous behavior leading sometime to new arc root position in back side and irregularities or disturbances on the arc root evolution versus time as observed at time t=18ms. In order to estimate the frequencies and speeds of rotation of the arc root and to compare them with the literature, we proceed as follows:

Knowing the number of revolutions 𝑁(𝑡) carried out by the arc root, given by the relation [START_REF] Cayla | Modélisation de l'interaction entre un arc électrique et une cathode[END_REF], we estimate instantaneous frequencies and rotational speeds of the arc root according to relations ( 29) and ( 30) respectively.

𝑁(𝑡) =

∑|𝜃(𝑡+∆t )-𝜃(𝑡)| 2𝜋 [START_REF] Cayla | Modélisation de l'interaction entre un arc électrique et une cathode[END_REF] Where:

-𝜃(𝑡) : represents the azimuthal position of the arc root at time 𝑡.

-𝜃(𝑡 + ∆t) : represents the azimuthal position of the arc root at time 𝑡 + ∆t -∆t =10 -5 s: is the calculation time step.

𝑓(𝑡) = N(t) 𝑛 𝑡 ∆t (29)
Where: 𝑛 𝑡 represents the number of time step ∆t associated to 𝑁(𝑡).

𝑣 𝑎𝑟𝑐 (𝑡) = 𝜋𝑑 𝑐 𝑓(𝑡)

𝑑 𝑐 is the internal cathode diameter. In the Figure 12, we present the frequency variations (left side) and the rotational arc root velocity versus time (right side), estimated respectively by relation ( 29) and [START_REF] Abramovitz | Handbook of mathematical functions[END_REF]. This figure shows that the frequency and the speed rotation of the arc increases rapidly from the first moments and then tends towards constant values of 1.14 kHz and 107m.s -1 respectively.

Figure 12 -Arc frequency and velocity versus time

Essiptchouk et al [START_REF] Essiptchouk | [END_REF] carried out an experimental study on the rotational speed of the electric arc within this torch technology. In their experiment, they deduce a rotational speed of the electric arc given by: 𝑣 𝑎𝑟𝑐 = 78 • 𝐼 4/9 𝐵 0.6 𝜌 0 -8/9 𝜑 -1/3 (31

)
Where 𝐼(𝐴), 𝐵(𝑇) and 𝜌 0 (kg/m 3 ) represent the current intensity, the magnetic field value and the initial mass density of the injected gas respectively. The parameter 𝜑 is given by the formulae:

𝜑 = 1 1+ 𝑣 𝑎 + 𝑣 𝑎 ( 32
)
𝑣 𝑎 is the gas velocity close to the arc. This study is performed for electric arc current intensities in air at atmospheric pressure between 100 and 1760A in the presence of an external magnetic field for magnetic field intensities between 0.005 and 3.9T, covering the power range of the torch studied in this paper.

Figure 13rotational arc root velocity versus axial gas velocity

From the equations ( 31) and [START_REF] Essiptchouk | [END_REF] given by Essiptchouk et al [START_REF] Essiptchouk | [END_REF] it is possible to represent the rotational arc root velocity versus the axial gas velocity. Several quantities are necessaries as the current intensity, the magnetic field, a mean mass density and the axial gas velocity. In order to draw the Figure 13, the current intensity is taken to 𝐼 =260A, the mass density as suggested by the author to its value at 𝑇 =300K (𝜌 0 =1.15 kg/m 3 ). An increase of the flow velocity leads to a decrease of the rotational velocity.

From the relation (31) of Essiptchouk et al, the axial velocity (𝑣 𝑎 ) is the most critical parameter to estimate. In their paper, the authors stipulate that the axial gas velocity is taken close to the arc but they haven't given any precision about its estimation. There for, in our model, we obtain for several values of inlet mass flow rate (D=18.5g/s, D=37g/s, D=74g/s) the axial gas velocity near the arc (see column n°6 of the Table 2). For the external magnetic field value, we have taken its maximum value at the plasma area interaction, B=200mT (Figures 9 and10). From the values of the axial gas velocity, the magnetic field, the current intensity and the mass density values, applying the equation [START_REF] Mutzke | IEEE Holm Conf. on Electrical Contacts[END_REF], it is possible to estimate a theoretical arc root rotational velocity by our model. In Table 2 and Figure 13, we present the arc root rotational velocity value obtained experimentally by Essiptchouk et al [START_REF] Essiptchouk | [END_REF] (see column n°7, in Table 2) and compared with the value founded by our model (last column in Table 2) estimated by relation [START_REF] Abramovitz | Handbook of mathematical functions[END_REF]. A good agreement is found between the rotational arc velocity found by the model and the formulation proposed by Essiptchouk. In this equation the torch dimensions as its radius are not explicit. Nevertheless, this dimension acts on the equation [START_REF] Essiptchouk | [END_REF], greater is the cathode radius, smaller the gas velocity and the parameter 𝜑 are. Another study from S. Ghorui et al [5] has shown that the arc velocity obeys to the relationship:

𝑣 𝑎𝑟𝑐  𝐵  (33)
In the literature the reported dependence of the coefficient varies from 0.34 to 0.60. In the case of an anode electrode and an arc current equals to 200A, with a magnetic field in reverse direction as us, the coefficient was found to 5]. In our study the same magnetic field value was applied for the three cases. Using and a value of the magnetic field equal to 200mT (2000 Gauss), the values of the arc velocities given in Table 2 lead from case 1 to 3 to proportionality coefficients respectively equal to 3.25, 2.65 and 2.58. For a given value of the magnetic field, the arc velocity increases with a diminution of the mass flow rate due to a decrease of the aerodynamic drag forces.

Nevertheless, it should be interesting as suggested by Brilhac et al [START_REF] Brilhac | [END_REF] to study if a correlation exists between the convection in the hollow cathode considering cold gas and the complete plasma flow representation in presence of an electric arc. From the initial injected mass flow rate D, one part goes down directly in the anode direction and the other part goes up in the hallow cathode before descending on the axis. The axial position where the integrated mass flow rate is 30% reduced corresponds to the beginning position of the arc attachment position found by the model. The study was performed for an inlet mass flow rate of D = 37 g/s and verified for two other values D = 18.5 g/s and D =74 g/s and a good agreement was found.

From these two empirics analysis the axial position and the rotational velocity can be deduced.

Plasma Characteristics

Based on the model developed, we present in this section general results on arc behavior. In the Figure 14, we present for different instants (10, 11 and 12ms), the temperature fields of the air plasma in middle plane (x=0) of © xxxx IOP Publishing Ltd the torch. At the anode exit (Torch exit, see Figure 1) we present in the Figures 15 and16, the fields and temperature profiles respectively for calculation instants, 10 to 15ms (step 1ms) in the two transverse directions (Ox) and (Oy). By observing the Figures 14, 15 and 16, we can notice that the fields and temperature profiles of the plasma show inhomogeneous distributions that evolve over time. Figure 17 shows the average temperature profiles in the two transverse directions (Ox) and (Oy) on the output plane (Torch exit, Figure 1) for different values of the integration time of 5, 10, 15, 20 and 24ms. We can see in the Figure 17 that the mean temperature profiles tend towards a Gaussian-type mean temperature profile with the integrated time starting from an integrated time of 15ms. Note that such a Gaussian temperature profile could be imposed in initial boundary conditions for models interested in the study of the variations of the thermal characteristics in the anode and at the torch exit.

Conclusion

This work was related to the study of a hollow cathode torch whose arc root is rotated under the influence of an external magnetic field. In order to apprehend the plasma characteristics, a 3D temporal model is developed in the upstream part of the torch describing the injectors, the flow in the cathode and in a plane located at the anode inlet. The magneto hydrodynamic equations are solved using the Ansys Fluent code. The model developed takes into account the magnetic field generated by the external coil. Thus the magnetic field resulting at each point is the contribution of the external field and the magnetic field generated by the flow of current in the electric arc. In order to describe the movement and the arc attachment position an effective conductivity module has been implemented. It allows the new position of the arc root to be determined temporally and to include an additional voltage drop without having to finely describe the non-equilibrium sheaths at the arcing point.

The model highlights the tortuous nature of the arc under the influence of convective forces and the resulting magnetic field. It clearly shows the existence of unsymmetrical plasma at the entrance of the downstream electrode. However if the fields are integrated over 15ms a Gaussian profile is observed.

© xxxx IOP Publishing Ltd

This tortuous character of the arc is also visible on the temporal variation of the tension. In relation to the parameters used, two characteristic times can be discerned: (1) the speed regulation time 𝑡 1 = 3ms, and (2) the stabilization time of the axial position, which is of the order of 𝑡 2 =26ms.

Our model does not take into account the erosion of the electrodes due to the interaction of the arc during its attachment or under the effect of radiation. However, it allows us to determine the axial position of the arc root and its rotation speed. In order to compare this attachment position with experimental work, Europlasma has carried out ultrasonic measurements to determine the wear zones relative to the attachment of the rotating arc root. A good correlation was observed between the plasma interaction area determined by the model and the erosion zone. In case of a cold flow study the arc attachment area corresponds to the position where the upstream flow was 30% reduced.

The rotation frequency determined by the model is in the order of 1.14 kHz after stabilization. The theoretical rotation speed is in the order of 102m/s. This theoretical velocity was compared to the one calculated from a formula in the literature. A good agreement was found validating the results of our model.

The model thus allows to characterize the plasma in the upstream electrode and to determine the spatial and temporal positions of the arc hooking. We will now use this model to describe the downstream part of the torch.

Figure 1 -

 1 Figure 1 -Overview of the studied device including the simplified scheme and position of the coil

Figure 2 -

 2 Figure 2 -Domains studied (Solid and Fluid) and grid

Figure 3 -

 3 Figure 3 -Initial conductor channel in the simulation

Figure 4 -

 4 Figure4-The layer defined for the Lindmayer model[17]

Figure 5 -

 5 Figure 5 -Field of the radial component of the magnetic field

Figure 6 -

 6 Figure 6 -Field of the axial component of the magnetic field Figures 5 and 6 show respectively the distributions of the radial and axial components of the magnetic field. The radial and axial magnetic field components have axisymmetric distributions. The maximum values of the radial and axial components of the magnetic field created by the coil are Brmax= 25mT and Bzmax=250mT, respectively. The black lines represent the position of the inner diameter of the cathode where the fluid flows.

Figure 7 -

 7 Figure 7 -Arc behavior within the geometry

Figure 8 -

 8 Figure 8 -Electric arc voltage versus time

Figure 9 -

 9 Figure 9 -Arc root axial position versus time in the left and, in the right radial (in blue) and axial (in red) external magnetic field components. I=260A, D=0.037kg/s.

  , we have represented the axial position of the © xxxx IOP Publishing Ltd

  rotation equilibrium zone predicted by the model, which we have compared to the experimental eroded zone obtained after several hours of running.

Figure 10 -

 10 Figure 10 -Wear position: experiment and model

Figure 11 -

 11 Figure 11 -Arc root azimuthal position versus time

Figure 14 -Figure 15 -Figure 16 -

 141516 Figure 14 -Temperature fields in the hallow cathode (t=10, 11 and 12ms)

Figure 17 -

 17 Figure 17 -Mean temperature profiles integrated on time at the exit plane for different integrated times on two perpendicular directions.

Table 2 -

 2 Confrontation of arc root rotational velocity between experiment[START_REF] Essiptchouk | [END_REF] and model

	Case	D(g/s)	I(A)	B(Gauss)	𝝆 𝟎 (kg/m 3 )	𝒗 𝒂 from the	𝒗 𝒂𝒓𝒄 [28]	𝒗 𝒂𝒓𝒄 from the
						model		model
	1	18.5	260	2000	1.15	20m/s	115 m/s	125 m/s
	2	37	260	2000	1.15	40 m/s	91.32 m/s	102 m/s
	3	74	260	2000	1.15	50 m/s	86 m/s	99.4 m/s

© xxxx IOP Publishing Ltd

Acknowledgement

The authors would like to thank the Europlasma Company for the ultrasonic measurements that were carried out.


Lucie Gallerand

Mathias Legrand

Thomas Dupont

Philippe Leclaire

Vibrational damping properties of finite microperforated plates

Microperforated plates (MPP) are traditionally used in sound control technologies to advantageously absorb acoustic waves. However, little is known concerning the additional structural damping they can induce through exchanges in the viscous and thermal boundary layers near the fluid-structure interface of the microperforations. MPP therefore offer an alternative to, or can be used together with, viscoelastic materials, commonly implemented to damp vibrations at medium and high frequencies. In this work, the structural damping capabilities of finite MPP are investigated. First, an analytical model is developed using an alternative form of the Biot model, classically devoted to porous plates, involving energy dissipation through viscous friction mechanisms. The analytical results are compared to experimental measurements of structural damping factors on various MPP samples. The model is validated and confirms the damping effect added by the microperforations in the low frequency range. A sensitivity analysis on perforation rate and perforation diameter provides a characteristic frequency that when coinciding with a plate natural frequency maximizes the added damping.

INTRODUCTION

Over the last few decades, microperforated plates (MPP) have been used as lightweight sound absorbers mainly because of their environment-friendly and hygienic properties. These simple structures can be made of various materials depending on the sought features: translucent (plexiglass), resistant to hostile environments (aluminum, steel, ...), biodegradable (wood, ...). These versatile sound absorption technologies are implemented, for instance, in cooling systems for turbines or in sound environments targeting comfort. These structures are known to dissipate acoustic energy as heat through exchanges in the viscous and thermal boundary layers near the fluid-solid interface of the microperforations. Multiple theoretical approaches have been developed to characterize the absorption coefficient of MPP. In this acoustic context, microperforated plates have been modelled by Maa [START_REF] Maa | Potential of microperforated panel absorber[END_REF] using Kirchhoff's equations. Maa's work has been extended in a vibroacoustic context where MPP acoustic properties have been also investigated in the light of structural vibration [START_REF] Takahashi | Flexural vibration of perforated plates and porous elastic materials under acoustic loading[END_REF][START_REF] Bravo | Vibroacoustic properties of thin micro-perforated panel absorbers[END_REF][START_REF] Lee | Widening the sound absorption bandwidths of flexible micro-perforated curved absorbers using structural and acoustic resonances[END_REF]. However, the contribution of the microperforated plate to vibration damping is generally neglected.

The structural dynamics of porous plates was concurrently explored using an alternative form of Biot's theory of poroelasticity [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range[END_REF][START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Higher frequency range[END_REF] leading to complicated equations. A simpler model in the form of two coupled partial differential equations governing the dynamics of a finite-size porous plate was proposed [START_REF] Leclaire | Transverse vibration of a thin rectangular porous plate saturated by a fluid[END_REF]. In this work, it showed that the porous plate induces additional damping at a characteristic frequency which is only function of porous plate geometrical parameters and fluid characteristics. Furthermore, an equivalent fluid concept already derived in for porous systems, was used in [START_REF] Atalla | Modeling of perforated plates and screens using rigid frame porous models[END_REF] for a sufficiently rigid structure to characterize acoustic proprertis of MPP. It has also been reported that fluid-solid interactions and visco-thermal effects induce additional damping in the case of MPP [START_REF] Lee | Sound absorption of a finite flexible micro-perforated panel backed by an air cavity[END_REF]. However, this effect has never been explored in detail.

The present work suggests to identify a finite size MPP to a porous plate with the aim of deriving an analytical model able to capture the additional damping features of MPP. In the first section, the analytical model is explained. The second section introduces the measurement procedure and corresponding results. Validation of the model is undertaken in a third section and closing comments are provided in the fourth and last section. The approaches and analyses proposed in the remainder are thoroughly developed in [START_REF] Gallerand | Vibration and damping analysis of a thin finite-size microperforated plate[END_REF], available in Open Access.

STRUCTURAL DYNAMIC OF FINITE-SIZE MPP

Analytical model

In this work, the microperforations have a circular cross-section of diameter d constant over the thickness h of the plate. By identifying the MPP, of dimensions L x × L y × h, to a porous plate, a model developed in the context of the vibration analysis of porous plates [START_REF] Leclaire | Transverse vibration of a thin rectangular porous plate saturated by a fluid[END_REF] is adapted to MPP. Bending only is considered and the MPP is regarded as a homogeneous plate excited by an external load f ext (x, y, t). The dynamic response of a MPP whose perforations are fluid-saturated is solution to the two coupled equations

D + α 2 M f h 3 12 ∇ 4 w s (x, y, t) + h(ρ ẅs (x, y, t) + ρ f ẅ(x, y, t)) = f ext (x, y, t), (1a)
αM f ∇ 2 w s (x, y, t) -ρ f ẅs (x, y, t) -m( jω) ẅ(x, y, t) = 0. ( 1b
)
Equation (1a) represents the elastic response of the homogeneous plate and Equation (1b) corresponds to the relative fluid-solid motion where

∇ 4 (•) = ∇ 2 (•) 2 and ∇ 2 (•) = ∂ 2 (•) ∂x 2 + ∂ 2 (•)
∂y 2 and w(x, y, t) = φ(w f (x, y, t) -w s (x, y, t)) with w f (x, y, t) is the fluid displacements and w s (x, y, t) the solid displacements. The parameters α and M f are elastic coefficients defined by Biot [START_REF] Biot | The elastic coefficients of the theory of consolidation[END_REF]. The coefficient D is the bending stiffness and the friction in the boundary layers in the perforations is represented by the complex coefficient

m( jω) = 1 φ ρ f α ∞ + σφ jω , (2)
where ρ f is the fluid density. The above definition of m( jω) is only effective valid in the low frequency range. Moreover, the parameters in the Biot theory (resistivity to airflow σ, perforation rate φ and tortuosity α ∞ ) are adapted to the particular case of MPP and listed in Table 1. Note that acoustic radiation in the perforation is considered in α ∞ through the parameter ε = 0.24 √ πd 2 (1 -1.14 √ φ) [START_REF] Atalla | Modeling of perforated plates and screens using rigid frame porous models[END_REF]. In Table 1, N perf is the number of microperforations and µ f is the dynamic viscoelasticity of the fluid. In addition, the mechanical parameters of the MPP (Young's modulus and solid mass density) Table 1: Geometrical parameter of a porous material adapted to MPP.

Airflow resistivity Perforated rate Tortuosity

σ = 32µ f φd 2 φ = N perf πd 2 4L x L y α ∞ = 1 + 2ε h
are also adapted to take account of perforations. Equivalent mechanical parameters are therefore defined [START_REF] Gallerand | Vibration and damping analysis of a thin finite-size microperforated plate[END_REF][START_REF] Boccaccini | A new approach for the Young's modulus-porosity correlation of ceramic materials[END_REF] as function of perforation diameter such that

E(φ) = E (1 -φ 2 ) 1 + (2 -3ν)φ , (3a)
ρ s (φ) = ρ s (1 -φ), (3b)
where E, ρ s and ν are respectively the Young's modulus, the solid density and the Poisson's coefficent of a non-perforated plate.

Characteristic frequency and maximum damping

The external load is set to zero and a modal analysis is performed to solve Equation (1) from which modal damping factors and the natural frequencies are readily obtained. An equivalent forced-response analysis is provided in [START_REF] Gallerand | Vibration and damping analysis of a thin finite-size microperforated plate[END_REF].

In the context of porous plate, a characteristic frequency f c is defined by Leclaire et al. [START_REF] Leclaire | Transverse vibration of a thin rectangular porous plate saturated by a fluid[END_REF] based on Biot's frequency. This frequency is characteristic of the viscous frictions between the fluid and the solid and depends only on the geometry of the porous materials and the parameters of the fluid. At f c , the additional damping induced by visco-thermic dissipation in the pores reaches a maximum [START_REF] Leclaire | Transverse vibration of a thin rectangular porous plate saturated by a fluid[END_REF]. By using MPP geometrical parameters defined in Table 1, the characteristic frequency is adapted to MPP as a function of d:

f c (d) = 32µ f 2πα ∞ ρ f d 2 . ( 4
)
In the following, MPP are made of aluminum and their loss factor is fixed at 0.01 %. The global loss factors therefore include both the structural damping of aluminum and the additional damping induced by the microperforations. Moreover, the natural resonance frequencies consider the structural damping and the damping added by the microperforations.

In Figure 1, the loss factor for the first MPP mode is plotted as a function of the dimensionless frequency f / f c (d) for three different values of φ. For a finite-size MPP, the additional damping is maximum for the mode mn when f c (d) coincides with a resonance frequency f mn where m and n are the indices of the modal discretization in the x and y directions. For a given φ, since the damping frequency is a function of d it is possible to adjust this last one to improve the damping capabilities of the MPP. Indeed, for a given mode there is a value of d for which the damping added by the microperforations is maximal. For a given φ, it appears in Figure 1 that the maximum damping is reached when there is coincidence between f 11 and f c (d). In other words, the additional damping provided by the microperforations is largest at the characteristic frequency, but also acts over a range of frequencies around it, as shown by the envelope in terms of damping around f / f c (d) = 1 in Figure 1. Moreover, for a given perforation diameter, the additional damping increases significantly with φ in Figure 1. The choice of φ is thus a compromise between the additional damping and the mechanical properties of the MPP since φ affects the natural frequencies by modifying the mass and the stiffness of the system.

In the following, the dimensions of the MPP along x and y are modified in order to explore the sensitivity of the first resonant frequency f 11 to the additional damping. Two different MPP, denoted MPP1* and MPP1, with constant φ are considered. For both MPP, d is fixed to induce a maximum damping on the first mode so that f c (d * comparison MPP1 * in Figure 2(a) and MPP1 in Figure 2(b) shows that the lower the natural frequency of the mode, the greater the damping added by the microperforations. This observation was confirmed numerically in many cases. The additional damping induced by the microperforations is therefore all the more important at low frequencies.

The influence of the considered vibration mode is now investigated. To this end, two MPP are considered. For the first one, denoted MPP1, the perforation diameter d ≡ d 1 is adjusted to induce maximum damping on the first vibration mode with natural frequency f c (d 1 ) = f 11 = 18 Hz. For the second one, denoted MPP2, the perforation diameter d ≡ d 2 is adjusted to induce maximum damping on the second vibration mode with natural frequency f c (d 2 ) = f 12 = 41 Hz. Their respective loss factor as a function of dimensionless frequency is illustrated in Figure 3. Both MPP share the same perforation rate φ = 10 % and plate dimensions. The comparison between Figure 3(a) and 3(b) shows that the maximum damping is observed for MPP1 and therefore for the first vibration mode.

With regard to Figure 2 and the damping will be all the more important as the first natural frequency f 11 is low.

EXPERIMENTAL VALIDATIONS

Experimental measurements are performed with an Oberst test bench illustrated in Figure 4. The MPP is clamped at x = 0, free at x = L x and free along y = ±L y /2. The excitation of amplitude F ext is provided by a magnetic exciter located at x = L x and the response is measured by a magnetic sensor at x = S . This set-up provides plate movements without adding extra mass and without obstructing the perforations, due to the sensors and exciters. The sample is designed with d = 1 mm to satisfy The mobility of the tested microperforated plate as well as the reference non-perforated plate is plotted in Figure 5. The decrease in the amplitude of the mobility of the MPP compared to that of the non-perforated plate shows the damping added by the microperforations, especially in the low frequency range. On the first peak, the gain is about 9 dB. The frequency shift is due to the φ-dependence of Young's modulus and solid mass density of the MPP [START_REF] Gallerand | Vibration and damping analysis of a thin finite-size microperforated plate[END_REF]. are probably due to the fact that Young's modulus is independent of the frequency in the analytical model.

CONCLUSIONS

The objective of this work was to investigate the possibility of using microperforated plates to archive structural damping. An analytical model used to explore the vibratory responses of MPP is developed from Biot's theory of poroelasticity. Results show that MPP could be used to significantly reduce structural vibrations in the low frequency range. This effect reaches a maximum at a particular frequency defined as a function of the perforation diameter and the perforation rate. The diameter of the perforation can be adjusted so that the characteristic frequency coincides with the natural frequency of the vibration mode to be damped. The additional damping induced by the microperforations is more important in the low frequency range. Therefore, it is recommended to adjust the perforation diameter to induce maximum damping on the first MPP mode.

5 %Figure 1 :

 51 Figure 1: Global loss factor for the first MPP mode for three different values of φ. Perforation diameter d set to maximize the additional damping on the first resonance, such that f c (d) = f 11 .

1 ) = f * 11 =

 111 10 Hz for MPP1 * and f c (d 1 ) = f 11 = 18 Hz for MPP1. The corresponding loss factor is plotted as a function of dimensionless frequency in Figure 2. The

Figure 2 :

 2 Figure 2: Global loss factor with constant φ = 10 %: (a) MPP size is 490 mm × 570 mm × 1 mm, (b) MPP size is 700 mm × 700 mm × 1 mm. In (a) and (b), d is set to induce maximum additional damping on the first resonance mode.

Figure 3 :

 3 Figure 3: Global loss factor as a function of dimensionless frequency for two MPPs with constant size (490 mm × 570 mm × 1 mm) and φ = 10 %: (a) first mode for MPP1, (b) second mode for MPP2. For both MPPs, d is set to induce maximum damping on the considered mode

Figure 4 :

 4 Figure 4: Experimental set-up: the dimensions are 130 mm × 30 mm × 0.87 mm; perforation rate and perforation diameter are d = 1 mm and φ = 10 % respectively.

Figure 5 :

 5 Figure 5: Mobilities for an MPP and the corresponding non-perforated plate.

lucie.gallerand.1@ens.etsmtl.ca