Nargisse Khiara

Fabien Onimus

Jean-Paul Crocombette

Laurent Dupuy

Thomas Pardoen

Jean-Pierre Raskin

Yves Bréchet

A molecular dynamics study of a cascade induced irradiation creep mechanism in pure copper

Keywords: Irradiation creep, Copper, Molecular Dynamics, Collision cascade, Shock wave Graphical

Recently, in-situ TEM straining experiments on pure copper have unraveled a high stress irradiation creep mechanism. Irradiation induced unpinning of dislocations from defects has been observed and the mean pinning lifetime has been determined. In the present study, molecular dynamics simulations are performed on pure copper to investigate the impact of collision cascades on screw dislocations pinned on Frank loops under high-applied stresses at 300 K in order to further quantitatively elucidate this mechanism. The simulations indicate two possible dislocation unpinning mechanisms. Unpinning can occur through loop destruction when the cascade is generated on the pinning points of the dislocation (type 1 unpinning). Unpinning can also be triggered by the shear stresses building up around a cascade generated in front of the dislocation in the glide plane (type 2 unpinning). Type 2 unpinning generally leads to dislocation repinning on cascade residues, so that it should only marginally contribute to irradiation creep. The mean pinning lifetime due to type 1 unpinning in the conditions of the in-situ TEM experiments is derived from a simple model previously developed for zirconium, and is found in the same orders of magnitude as in the experiments.

Introduction

Irradiation creep deformation is a visco-plastic deformation phenomenon occurring in components inside nuclear reactors at temperatures below 0.3Tm under applied load and irradiation. Irradiation creep has been mainly investigated through in-reactor experiments, allowing for the measurement of creep parameters and their dependence on neutron flux and fluence, temperature, stress, and material characteristics [START_REF] Aitkhozhin | Radiation-induced creep of copper, aluminium and their alloys[END_REF][START_REF] Garnier | Irradiation creep of SA 304L and CW 316 stainless steels: Mechanical behaviour and microstructural aspects. Part I: Experimental results[END_REF][START_REF] Onimus | Irradiation Creep in Materials[END_REF][START_REF] Singh | Final report on in-reactor tensile tests on OFHC-Copper and CuCrZr alloys[END_REF]. However, such experiments are not amenable to the in-situ characterization of the elementary physical mechanisms at the origin of the macroscopic behavior. Therefore, although the macroscopic phenomenology for irradiation creep is relatively well rationalized, the underlying microscopic mechanisms are still unclear [START_REF] Matthews | Irradiation creep models-An overview[END_REF][START_REF] Garner | Irradiation creep mechanisms: an experimental perspective[END_REF]. Many theoretical mechanisms have been proposed in the literature, all involving atomic displacements and generation of point defects under irradiation, but differing in the detailed processes of their redistribution under applied stress [START_REF] Onimus | Irradiation Creep in Materials[END_REF][START_REF] Matthews | Irradiation creep models-An overview[END_REF][START_REF] Wolfer | Correlation of radiation creep theory with experimental evidence[END_REF][START_REF] Mansur | Irradiation creep by dislocation glide enabled by preferred absorption of point defects-theory and experiment[END_REF].

Very recently, a high stress irradiation creep phenomenon, first unraveled for a zirconium alloy was extended to pure copper which is considered as a model material for face-centered cubic (FCC) structures. Through in-situ transmission electron microscopy straining experiments under heavy ion irradiation, the authors [START_REF] Gaumé | In-situ observation of an irradiation creep deformation mechanism in zirconium alloys[END_REF][START_REF] Khiara | Insitu TEM irradiation creep experiment revealing radiation induced dislocation glide in pure copper[END_REF] observed an irradiation induced unpinning of dislocations from irradiation defects followed by dislocation glide. They proposed that the dislocation unpinning from irradiation defects is triggered by displacement cascades generated by irradiation. The mean pinning lifetime, i.e., the time necessary for a dislocation to unpin from an irradiation defect, was estimated around 25 and 40 seconds at room temperature for the zirconium alloy and pure copper, respectively. Molecular dynamics (MD) is instrumental to investigate elementary mechanisms at the atomic level. MD simulations in copper related to the present context focused on dislocation interaction with defects [START_REF] Doyama | Interactions between point defects and edge dislocations in copper[END_REF][START_REF] Lee | Molecular dynamics simulation of screw dislocation interaction with stacking fault tetrahedron in face-centered cubic Cu[END_REF][START_REF] Nogaret | Atomic-scale plasticity in the presence of Frank loops[END_REF][START_REF] Osetsky | Atomic-scale study of dislocation-stacking fault tetrahedron interactions. Part I: mechanisms[END_REF][START_REF] Rodney | Molecular dynamics simulation of screw dislocations interacting with interstitial frank loops in a model FCC crystal[END_REF] and cascade damage [START_REF] Calder | On the origin of large interstitial clusters in displacement cascades[END_REF][START_REF] De La Rubia | Molecular dynamics studies of defect production and clustering in energetic displacement cascades in copper[END_REF]. The effect of displacement cascades on straight dislocations was studied for aluminum [START_REF] Voskoboinikov | Interaction of collision cascades with an isolated edge dislocation in aluminium[END_REF][START_REF] Voskoboinikov | MD simulations of collision cascades in the vicinity of a screw dislocation in aluminium[END_REF], tungsten [START_REF] Fu | Effect of collision cascades on dislocations in tungsten: A molecular dynamics study[END_REF], Fe-20Cr-25Ni alloy [START_REF] Lu | Interactions between displacement cascade and dislocation and their influences on Peierls stress in Fe-20Cr-25Ni alloys[END_REF], and Fe [START_REF] Heredia-Avalos | Collision cascade effects near an edge dislocation dipole in alpha-Fe: Induced dislocation mobility and enhanced defect clustering[END_REF]. Some studies [START_REF] Voskoboinikov | Interaction of collision cascades with an isolated edge dislocation in aluminium[END_REF][START_REF] Fu | Effect of collision cascades on dislocations in tungsten: A molecular dynamics study[END_REF][START_REF] Lu | Interactions between displacement cascade and dislocation and their influences on Peierls stress in Fe-20Cr-25Ni alloys[END_REF] evidenced that, when a cascade is generated on an edge dislocation, the dislocation can capture the generated point defects, leading to dislocation climb. For a screw dislocation, the dislocation can either capture point defects and form a helical turn, or, if subjected to loading, cross-slip [START_REF] Voskoboinikov | MD simulations of collision cascades in the vicinity of a screw dislocation in aluminium[END_REF][START_REF] Lu | Interactions between displacement cascade and dislocation and their influences on Peierls stress in Fe-20Cr-25Ni alloys[END_REF]. An intriguing cascade effect on dislocations was first observed by Fu et al. [START_REF] Fu | Effect of collision cascades on dislocations in tungsten: A molecular dynamics study[END_REF] in tungsten: when a cascade is generated near an edge dislocation, kinks are formed near the center of the cascade, and the dislocation glides toward the cascade. A similar phenomenon was observed by Heredia-Avalos et al. [START_REF] Heredia-Avalos | Collision cascade effects near an edge dislocation dipole in alpha-Fe: Induced dislocation mobility and enhanced defect clustering[END_REF] in iron. The authors observed that a cascade generated near an edge dislocation dipole leads to the glide of the dislocations toward the cascade center. The authors reproduced the same phenomenon by generating a thermal spike near the dislocation dipole. Heredia-Avalos et al. [START_REF] Heredia-Avalos | Collision cascade effects near an edge dislocation dipole in alpha-Fe: Induced dislocation mobility and enhanced defect clustering[END_REF] stated that the observed dislocation motion is triggered by the shock wave generated by the cascade (or the thermal spike). Starostenkov et al. [START_REF] Starostenkov | Dynamics of edge dislocations in a low-stability FCC-system irradiated by high-energy particles[END_REF] investigated the post cascade shock wave effect on a single edge dislocation or on an ensemble of edge dislocations in nickel, and the authors observed the glide of the dislocations toward the wave source. Korchuganov et al. [START_REF] Korchuganov | Mobility of edge dislocations in stressed iron crystals during irradiation[END_REF][START_REF] Korchuganov | MD simulation of primary radiation damage in metals with internal structure[END_REF] found that a shock wave generated near an edge dislocation under applied stress in iron induces dislocation glide toward the wave source. The authors concluded that this phenomenon can contribute to irradiation creep.

In a previous study, we [START_REF] Khiara | A novel displacement cascade driven irradiation creep mechanism in α-zirconium: A molecular dynamics study[END_REF] investigated the effect of a cascade on a screw dislocation pinned on a loop at high applied stress in zirconium using MD simulations to reproduce similar configurations as in the TEM in-situ straining experiments of Gaumé et al. [START_REF] Gaumé | In-situ observation of an irradiation creep deformation mechanism in zirconium alloys[END_REF] on a zirconium alloy. A possible unpinning of dislocations by cascades was predicted at high stress levels slightly inferior to the critical stress for unpinning (without a cascade). In these simulations, the generated cascades had to encompass the pinning points of dislocations to induce unpinning, and we developed a simple closed-form model to account for this phenomenon. Based on this closed-form model, we [START_REF] Khiara | A novel displacement cascade driven irradiation creep mechanism in α-zirconium: A molecular dynamics study[END_REF] estimated the mean pinning lifetime in the conditions of the experiments of Gaumé et al. [START_REF] Gaumé | In-situ observation of an irradiation creep deformation mechanism in zirconium alloys[END_REF] and obtained pinning lifetimes at high stress levels in the same orders of magnitude as the experimental pinning lifetime.

The objective of this work is to extend the numerical study using MD simulations made for zirconium to pure copper and to compare the MD results to our recently published experimental results on pure copper [START_REF] Khiara | Insitu TEM irradiation creep experiment revealing radiation induced dislocation glide in pure copper[END_REF].

MD simulation method 2.1. Interaction dislocation / loop

All the molecular dynamics simulations are performed on pure copper (FCC structure) using the Largescale Atomic and Molecular Massively Parallel Simulator (LAMMPS) [START_REF] Plimpton | Fast Parallel Algorithms for Short-Range Molecular Dynamics[END_REF]. The embedded atom method (EAM) used is a potential derived and updated by Ackland et al. [START_REF] Ackland | Simple N-body potentials for the noble metals and nickel[END_REF] to take into account close-range repulsion for irradiation studies (Cu1_v2). Interactions between screw dislocations and irradiation defects were observed during our previous in-situ straining experiments [START_REF] Khiara | Insitu TEM irradiation creep experiment revealing radiation induced dislocation glide in pure copper[END_REF]. The objective is to reproduce a similar configuration as in the experiments. The , , and axis of the simulation box correspond respectively to the directions 110 , 1 12 and 11 1 . The box dimensions are approximatively 400 Å × 234 Å × 220 Å, along the , and directions respectively, which corresponds to around 1.6 million atoms. The simulation box dimension along the axis is close to the mean obstacle spacing in our in-situ straining experiments under irradiation [START_REF] Khiara | Insitu TEM irradiation creep experiment revealing radiation induced dislocation glide in pure copper[END_REF] (estimated between 45 and 70 nm). Periodic boundaries are applied in the and directions.

The methodology is similar to the one used for zirconium in our previous MD study [START_REF] Khiara | A novel displacement cascade driven irradiation creep mechanism in α-zirconium: A molecular dynamics study[END_REF]. First, a molecular dynamics simulation box containing a screw dislocation and a Frank loop is generated using NUMODIS [START_REF] Drouet | A direct comparison between in-situ transmission electron microscopy observations and Dislocation Dynamics simulations of interaction between dislocation and irradiation induced loop in a zirconium alloy[END_REF][START_REF] Shi | Interaction of< 1 0 0> dislocation loops with dislocations studied by dislocation dynamics in α-iron[END_REF][START_REF]PROJECT NUMODIS : Numerical Modelling of Dislocations[END_REF]. NUMODIS is a Dislocation Dynamics simulation code, but has recently incorporated a tool to generate simultaneously MD and DD identical simulation boxes to perform similar MD and DD simulations.

The screw dislocation has a Burgers vector = √ 110], and a glide plane (11 1). The Frank loop has a (1 1 1 ) habit plane and contains 61 SIAS (self-interstitial atoms), which corresponds to an approximate loop diameter of 2.3 nm. The loop dimension properly compares to the 2 nm loops observed in our experiments [START_REF] Khiara | Insitu TEM irradiation creep experiment revealing radiation induced dislocation glide in pure copper[END_REF]. The Frank loop center lies in the dislocation glide plane. As for the previous MD study on zirconium, this specific configuration is chosen since the interaction of the chosen dislocation and Frank loop leads to the formation of a helical turn (Figure 1.c), which is the configuration for which the obstacle is the strongest [START_REF] Nogaret | Atomic-scale plasticity in the presence of Frank loops[END_REF].

The simulation box is first thermally equilibrated at 300 K for 50 ps with a time-step of 5 fs. The crystal is then divided into three layers along the Z direction. The atoms of the top and bottom layers are held in fixed positions with respect to one another, while no constraint is imposed on the atoms of the middle layer. As presented in Figure 1.a.1, a shear strain rate of 10 ps -1 is applied on the atoms of the top layer (in green), while the atoms of the bottom layer are held in fixed positions (in grey). The corresponding shear stress is derived from the force applied by the free atoms of the middle layer on the upper layer atoms. The shear stress evolution with strain is presented in Figure 1.b. In order to smooth out stress fluctuations (blue curve), an averaging over 1000 time-steps is performed (red curve), corresponding to a duration of 250 ps and a strain increment of 0.025%. The critical stress for dislocation unpinning is the shear stress necessary for the release of the dislocation (without cascade) and is equal to 200 MPa.

The processing and visualization tool OVITO [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with OVITO-the Open Visualization Tool[END_REF] is used to analyze and post-process the MD results. The common neighbor analysis tool (CNA) [START_REF] Honeycutt | Molecular dynamics study of melting and freezing of small Lennard-Jones clusters[END_REF] is used to display atoms that are not in perfect FCC positions, whereas the dislocation extraction algorithm (DXA) [START_REF] Stukowski | Extracting dislocations and non-dislocation crystal defects from atomistic simulation data[END_REF] is used to display dislocations. The Voronoi analysis method is used to calculate the atomic volume with OVITO [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with OVITO-the Open Visualization Tool[END_REF][START_REF] Béland | Atomistic material behavior at extreme pressures[END_REF]. This allows calculating the atomic stresses in units of pressure from LAMMPS atomic stress outputs.

Cascade effect on pinned dislocation

In order to investigate the effect of a cascade on a screw dislocation pinned on a Frank loop for stresses slightly inferior to the critical stress for unpinning, the same methodology as the one proposed for zirconium is used [START_REF] Khiara | A novel displacement cascade driven irradiation creep mechanism in α-zirconium: A molecular dynamics study[END_REF]. Four high stress levels below the critical stress are selected: 0.85 , 0.9 , 0.95 and 0.975 (white points in Figure 1.b). The atomic configurations corresponding to the selected shear stresses are extracted (Figure 1.a.2). After a thermal equilibration at 300 K for 50 ps with a time-step of 5 fs, the atoms of the top and bottom layers are held fixed (in grey in Figure 1.a.2), so that the total strain is kept constant. To generate a displacement cascade, an energy (between 10 and 30 keV) along the direction is transferred to an atom, called Primary Knock-on-Atom (PKA). Since the PKA velocity direction is along the direction, the center of the cascade generated lies over the PKA initial position along the axis. A variable time step is defined in order to limit the distance traveled by atoms to less than 1% of the lattice constant. The simulations are run for 150 000 timesteps, corresponding to approximately 100 ps, which is long enough to ensure that ballistic and thermally enhanced recovery stages of the cascade are over. The temperature is set at 300 K in order to allow for comparisons with our in-situ experiments performed at room temperature [START_REF] Khiara | Insitu TEM irradiation creep experiment revealing radiation induced dislocation glide in pure copper[END_REF].

In this study, the impact of the shear stress level, PKA position, and energy are investigated. Since various configurations are studied, the PKA positions and energies selected for each configuration will be detailed in the following sections.

Results

Impact of the PKA position

The unpinning of screw dislocations from Frank loops is predicted in the two different configurations presented in Figure 2. The first configuration, called type 1 unpinning, is similar to the one we observed for zirconium [START_REF] Khiara | A novel displacement cascade driven irradiation creep mechanism in α-zirconium: A molecular dynamics study[END_REF]. In this case, the cascade encompasses the helical turn. After the end of the cascade, the dislocation is free and the loop is destroyed. Residual clusters are formed behind the dislocation. Unlike in zirconium for which the loop is only partially destroyed, the Frank loop in copper is completely destroyed [START_REF] Khiara | A novel displacement cascade driven irradiation creep mechanism in α-zirconium: A molecular dynamics study[END_REF].

The second configuration that induces dislocation unpinning from the loop, called type 2 unpinning, occurs for cascades generated in front of the dislocation along the glide direction. As shown in Figure 2, when the cascade peak occurs, the cascade starts to attract one or both segments of the dislocation.

The dislocation segments enter in contact with the displaced atoms within the cascade. The dislocation then unpins from the initial Frank loop, leaving the Frank loop unfaulted, as in the case of a standard dislocation unpinning from a Frank loop with increasing stress without cascade [START_REF] Nogaret | Atomic-scale plasticity in the presence of Frank loops[END_REF]. The dislocation is often repinned on residual clusters left after cascade relaxation and forms a helical turn (Figure 2), as observed by Voskoboinikov et al. [START_REF] Voskoboinikov | MD simulations of collision cascades in the vicinity of a screw dislocation in aluminium[END_REF] for displacement cascades generated near straight screw dislocations in aluminum. Type 2 unpinning is thus divided into two subcategories: unpinningrepinning mechanism (when dislocations unpin from loops but repin on cascade residues), and permanent unpinning (when dislocations remain straight and free to glide under stress). Permanent unpinning generally occurs when the cascade is generated slightly below the dislocation glide plane (Figure 2).

To ensure that the two unpinning types do not result from box size effects, calculations were repeated for a larger simulation box size along the axis (405 Å rather than 234 Å), and both unpinning types, with similar trends, were captured for the larger box in comparison with the smaller box. More details on this analysis can be found in Supplementary Data. Likewise, we checked using various PKA directions that both unpinning can occur when PKA initial velocity is along axes other than the axis.

Impact of stress

The impact of the stress level on the unpinning probability is investigated at different PKA positions. As stated earlier, four stress levels are studied: 0.85 , 0.9 , 0.95 and 0.975 . A PKA position grid around the dislocation is addressed. For 0.85 and 0.9 , the grid is of 160 Å × 160 Å × 160 Å, with a spacing between PKA positions of 20 Å. For 0.95 and 0.975 , a grid of 160 Å × 220 Å × 160 Å is selected with a spacing between PKAs of 20 Å. The PKA energy is set at 20 keV and the PKA direction is always along the axis.

The grids and simulation results are presented in Figure 3: the small blue points represent the PKA positions that did not induce unpinning, the green and red points represent PKA that triggered type 1 and type 2 unpinning, respectively. Increasing unpinning frequency with increasing stress is observed. Type 1 unpinning prevails at lower stresses (0.85 and 0.9 , ) and type 2 unpinning prevail at higher stresses (0.95 and 0.975 ). This dependency on the sensitivity to stress is more important for type 2 unpinning. Moreover, the repinning probability is around 94%("4%) at 0.95 and around 80%("4%) at 0.975 . Therefore, it appears that type 2 permanent unpinning probability increases with stress.

To better visualize the PKA position effect on unpinning, the unpinning probability at both 0.95 and at 0.975 projected along the , and directions are represented in Figure 4 and Figure 5, respectively. For type 1 unpinning, over the and directions, the maximum unpinning probability is centered on the helical turn. For type 2 unpinning, over the direction, the unpinning probability is at its maximum around 50 # 60 Å in front the mean dislocation position along the axis. Over the direction, the maximum probability is around the helical turn, but with a wider spread in comparison with type 1 unpinning. The wider spread for type 2 unpinning along the axis can be explained by the fact that the cascade can trigger dislocation unpinning by only attracting one of the segments of the dislocation. Over the axis, the maximum unpinning probability is observed approximatively at the same position for both unpinning types. The maximum is observed under the dislocation along the axis, since the PKAs at these positions generate cascades approximately on the dislocation glide plane.

To summarize, type 1 unpinning occurs when the cascade is on the helical turn, whereas type 2 unpinning takes place when the cascade is generated in front of the dislocation in the glide plane.

Impact of PKA energy

The impact of the PKA energy is investigated at a selected stress level of 0.95 . Four PKA energies are investigated: 10 keV, 15 keV, 20 keV, and 30 keV. The PKA position grid investigated is of 160 Å × 160 Å × 160 Å with a spacing between PKA positions of 40 Å. The PKA velocity direction is always along the axis. The PKA position grid and simulation results are presented in Figure 6. Increasing unpinning probability is observed with increasing PKA energy. The sensitivity to the PKA energy seems higher for type 2 unpinning. The repinning probability is equal to 100% at 10 keV, 83% at 15 keV, 86% at 20 keV and 80% at 30 keV. There are, however, large uncertainties on these values since there is only 1 effective type 2 unpinning at 10 keV, 6 at 15 keV, 7 at 20 keV, and 10 at 30 keV. Therefore, more investigation is needed to draw any conclusions on the effect of the PKA energy on the repining probability.

Discussion 4.1. Type 2 unpinning

Type 2 unpinning occurs for cascades generated in front of the dislocation along the glide direction through a process similar to an increase in stress, with an unfaulted loop left after the unpinning process. When the cascade peak occurs, one or both segments of the dislocation glide toward the cascade center and enter in contact with the displaced atoms within the cascade. The dislocation then unpins from the Frank loop.

Type 2 unpinning resembles behavior seen in previous work on cascade, thermal spike, and shock wave effects on straight dislocations. The dislocation motion triggered by a cascade (or a shock wave) toward the cascade center (or the wave source) was reported by several authors for different metals through MD simulations [START_REF] Fu | Effect of collision cascades on dislocations in tungsten: A molecular dynamics study[END_REF][START_REF] Heredia-Avalos | Collision cascade effects near an edge dislocation dipole in alpha-Fe: Induced dislocation mobility and enhanced defect clustering[END_REF][START_REF] Korchuganov | MD simulation of primary radiation damage in metals with internal structure[END_REF][START_REF] Korchuganov | Mobility of edge dislocations in stressed iron crystals during irradiation[END_REF][START_REF] Starostenkov | Dynamics of edge dislocations in a low-stability FCC-system irradiated by high-energy particles[END_REF]. Cascade induced glide of dislocations was also theoretically investigated by Trinkaus [START_REF] Trinkaus | Local stress relaxation in thermal spikes as a possible cause for creep and macroscopic stress relaxation of amorphous solids under irradiation[END_REF] and Zhukov and Boldin [START_REF] Zhukov | Interaction of Cascade-Generated Nonlinear Elastic Waves with Dislocation Structure Elements in Metals[END_REF].

Heredia-Avalos [START_REF] Heredia-Avalos | Collision cascade effects near an edge dislocation dipole in alpha-Fe: Induced dislocation mobility and enhanced defect clustering[END_REF] attributed the initial dislocation motion toward the displacement cascade to the shock wave generated by the cascade. The same observation is made in our MD simulations. Figure 7 shows frames representing the initial dislocation motion during type 2 unpinning, along with the atomic pressure and shear stress fields. A correlation is observed between the dislocation motion and the time of arrival of the shear stress shock wave (and the pressure shock wave) at the dislocation. As reported by Béland et al. [START_REF] Béland | Atomistic material behavior at extreme pressures[END_REF] , the shock wave appears during the sonic step of the cascade, and propagates during the sonic and thermal annealing steps. As shown in Figure 7, the dislocation starts its first motion when the shock wave arrives near the dislocation. There is also a correlation between the first dislocation motion and the sign of the shear stress wave. At a time of about 1.21 ps and 1.78 ps in Figure 7, the right segment moves significantly toward the cascade center whereas the left segment moves slightly backward. This can be simply linked to the signs of the shear stresses and to the resultant Peach-Köhler force.

Far from the dislocation, the same shear stress distribution as in perfect crystal is observed: negative shear stress in the upper right and positive shear stress in the upper left (Figure 7). Near the dislocation, the negative shear stress wave is either attenuated or absent (Figure 7) in comparison with a perfect crystal (see supplementary data for more examples). Therefore, in the area of negative shear stress, the dislocation can either partially and slightly move backward (Figure 7) or is unaffected. On the contrary, in the positive shear stress area, the dislocation is highly affected by the stress shock wave and partially glides toward the center of the cascade.

The maximum unpinning probability is found for initial PKA positions at about 50 # 60 Å from the mean dislocation position along the axis. This can easily be explained in terms of cascade radius (Figure 4 and Figure 5). Indeed, the shear stresses are generated around the cascade core, which, at its peak, has a radius of 40 # 50 Å for 20 keV cascades.

Type 2 unpinning probability increases with increasing PKA energy and increasing stress. This is analogous to the observations of Korchuganov et al. [START_REF] Korchuganov | Mobility of edge dislocations in stressed iron crystals during irradiation[END_REF] for straight edge dislocations under applied shear stress and exposed to a shock wave in iron. The authors observed that the larger the shear stress and the amplitude of the shock wave, the larger the displacement of the dislocation.

One of the main objectives of this work is to estimate, based on the MD calculations, the mean pinning lifetime in the experimental conditions of Ref. [START_REF] Khiara | Insitu TEM irradiation creep experiment revealing radiation induced dislocation glide in pure copper[END_REF]. After dislocation unpinning from the Frank loop, there are two possible outcomes: the dislocation is either repinned (unpinning-repinning process) or completely free (permanent unpinning process). Type 2 unpinning appears and becomes predominant at very high stress level. However, there is a very high repinning probability on cascade residues. For 20 keV cascades, up to a stress level of 0.95 , permanent type 2 unpinning represents less than 7% of type 1 unpinning and is therefore negligible.

For 0.975 , type 2 permanent unpinning induces an increase in the unpinning probability of nearly 50%. Therefore, except for very high stress levels, type 2 permanent unpinning is negligible in comparison with type 1 unpinning. Moreover, in Ref. [START_REF] Khiara | Insitu TEM irradiation creep experiment revealing radiation induced dislocation glide in pure copper[END_REF], the dislocation mean free path after unpinning is estimated at around 45 nm, which is not compatible with type 2 unpinning-repinning mechanism. Indeed, the mean dislocation glide distance before repinning is around 0.7 and 0.6 nm with a maximum traveled distance of 7 and 8 nm at respectively 0.95 and 0.975 , which is way below the experimental dislocation free path.

Another possible effect of type 2 unpinning on irradiation creep can arise from the potential lower obstacle strength of cascade residues in comparison with the initial Frank loop. In order to investigate this effect, 25 post-cascade MD straining calculations are performed on configurations for which type 2 unpinning-repinning process was triggered by a 20 keV cascade. The post unpinning-repinning critical stress is found at around 93% of the initial critical stress. More details on this analysis can be found in supplementary data. This shows that type 2 unpinning-repinning might be equivalent to a very slight reduction in obstacle strength. However, there is a wide spread of critical stress over and under the initial critical stress. Therefore, more investigations on this subject, especially on the effect of loop size and PKA energy, are needed to draw any definitive conclusion.

Type 1 unpinning

The type 1 unpinning process is driven by the destruction of the loop. Following the conclusions of the previous analysis on type 2 unpinning, type 1 unpinning has a predominant effect on irradiation creep, at least up to 0.95 .

Based on the MD calculations, the unpinning probability and standard deviation are defined as follows:

+ = , -./0 , -1- 2 , ( 1
)
3$ = 4 +(1 # +) , -1- 5 . (2)
where , -./0 is the number of effective type 1 PKAs and , -1-is the number of investigated points.

From the unpinning probability, an effective volume and the corresponding standard deviation can be defined through:

* 066 = + × * -1-, (3)
3$ 7 899 = 3$ × * -1-. (4)
As stated earlier, type 1 unpinning is triggered when the cascade encompasses the helical turn. In order to take this observation into account, we proposed a simple model to calculate the effective volume for zirconium: if a cascade with a radius : , that depends on the PKA energy, encompasses the pinning points of the dislocation, type 1 dislocation unpinning occurs. The distance between pinning points ∆ <11/ depends on stress. The higher the stress, the lower this distance.

Therefore, the effective volume can be simply calculated as [START_REF] Khiara | A novel displacement cascade driven irradiation creep mechanism in α-zirconium: A molecular dynamics study[END_REF]:

if : ($ %&' ) < ∆ >??@ (A) * 066 ($ %&' , ) = 0 , ( 5
) if : ($ %&' ) > ∆ >??@ (A) * 066 ($ %&' , ) = D E4: ($ %&' ) ∆ <11/ ( )F E2: ($ %&' ) # ∆ <11/ ( )F . (6)
As in Ref. [START_REF] Khiara | A novel displacement cascade driven irradiation creep mechanism in α-zirconium: A molecular dynamics study[END_REF], to calculate cascade radius for a given PKA energy, 10000 cascade calculations are performed using Iradina [START_REF] Borschel | Ion beam irradiation of nanostructures -A 3D Monte Carlo simulation code[END_REF][START_REF] Crocombette | Quick calculation of damage for ion irradiation: implementation in Iradina and comparisons to SRIM[END_REF] for the investigated PKA energy (here 20 keV). Moreover, the pinning point spacing positions is taken to be the distance between the upper and the bottom parts of the loop (as shown by the insert in Figure 8.a).

Figure 8.a represents the evolution of the effective volume with stress for 20 keV cascades based on the MD calculations (in green) and on the simple model developed in Ref. [START_REF] Khiara | A novel displacement cascade driven irradiation creep mechanism in α-zirconium: A molecular dynamics study[END_REF] (in red). The error bars for the model come from the dispersion in cascade radii. The model seems to be in good agreement with MD calculations, although a slight underestimation is observed.

Based on the hypothesis that the effective volume around dislocation pinning points for a given stress and PKA energy can be derived from Eq. 5 and 6, we [START_REF] Khiara | A novel displacement cascade driven irradiation creep mechanism in α-zirconium: A molecular dynamics study[END_REF] developed a closed-form model to estimate the mean pinning lifetime evolution with stress. This model accounts for cascade fragmentation into subcascades. For copper, the threshold energy for cascade fragmentation into subcascades is found to be around 60 keV [START_REF] Backer | A model of defect cluster creation in fragmented cascades in metals based on morphological analysis[END_REF].

Following the model developed before, the evolution of the pinning lifetime with stress in the experimental conditions of the in-situ straining experiments under irradiation in Ref. [START_REF] Khiara | Insitu TEM irradiation creep experiment revealing radiation induced dislocation glide in pure copper[END_REF] is represented in Figure 8.b. At high stress levels, the pinning lifetimes are in the same ranges as the experimental value we obtained [START_REF] Khiara | Insitu TEM irradiation creep experiment revealing radiation induced dislocation glide in pure copper[END_REF]. Therefore, the proposed cascade induced unpinning mechanism can provide an explanation to the irradiation triggered unpinning of dislocations at high stress levels (over 0.8 ) observed experimentally. Moreover, the minimum stress for cascade induced unpinning can also be derived from this model. This minimum stress is found equal to 0.76 and corresponds to the stress level for which the condition : ($ %&' ) < ∆ >??@ (A) is valid for the largest possible cascades (60 keV cascades given the threshold energy for cascade fragmentation).

Conclusion

MD simulations are performed on pure copper, considered as a model material for FCC metals, in order to investigate the possible unpinning by a displacement cascade of a screw dislocation pinned on a Frank loop and provide new insights on high stress irradiation creep mechanisms.

The main findings of this work are the following.

• The possible unpinning of screw dislocations pinned on a Frank loop occurs at high stress levels (at least starting from 0.85 ); • Two unpinning processes are observed: type 1 and type 2 unpinning; • Type 1 unpinning is similar to the one we have previously unraveled for zirconium [START_REF] Khiara | A novel displacement cascade driven irradiation creep mechanism in α-zirconium: A molecular dynamics study[END_REF]: the cascade occurs on the helical turn and induces the destruction of the loop, leaving a free dislocation. • Type 2 unpinning is triggered by the shear stresses generated around the cascade core. The dislocation unpins from the loop as in the case of an increase in shear stress: the loop is left unfaulted. In a large majority of cases, the dislocation is very often pinned on cascade residues (unpinning-repinning process). In very few cases, the dislocation is free (permanent unpinning process); • Type 2 unpinning occurs and becomes dominant at very high stress levels (> 0.95 ); • Both unpinning types are more likely to occur with increasing stress and increasing PKA energy;

• Since type 2 unpinning is observed at very high stress levels and leads in the majority of cases to dislocation repinning, its impact on irradiation creep is most likely far smaller than type 1 unpinning; • The model we developed previously for zirconium [START_REF] Khiara | A novel displacement cascade driven irradiation creep mechanism in α-zirconium: A molecular dynamics study[END_REF] is applied to reproduce type 1 unpinning: if the cascade encompasses the helical turn, the dislocation unpins. This model slightly underestimates the unpinning probability but gives values in the same orders as the ones obtained from the MD calculations. Using this simple model, the mean pinning lifetime of dislocations is derived and values similar to the experimental pinning lifetime estimated in Ref. [START_REF] Khiara | Insitu TEM irradiation creep experiment revealing radiation induced dislocation glide in pure copper[END_REF] are obtained. Therefore, this cascade induced unpinning of dislocation process may provide an explanation for the high stress irradiation creep phenomenon observed during insitu experiments.

Acknowledgment

This work was granted access to the HPC resources of TGCC under the allocations 2019-A0070911020 and 2020-A0090911020 made by GENCI.

Figure 1 -

 1 Figure 1 -Boundary conditions (using the CNA method): (a.1) for the interaction of the dislocation with the loops -(a.2) for the cascade calculations -(b) Stress-strain curve for the interaction dislocation loop (without cascade) before (blue) and after (red) curve smoothing, the white points represent the shear stresses investigated: 0.85 , 0.9 , 0.95 and 0.975 -(c) Helical turn displayed using the DXA algorithm.

Figure 2 -

 2 Figure 2 -Snapshots showing the time evolution of core atoms (using the CNA method) in a case where type 1 and type 2 unpinning processes of the dislocation are triggered for 20 keV cascades (the corresponding videos can be found in the Supplementary data).

Figure 3 -

 3 Figure 3 -MD simulation results for four stress levels for a PKA energy $ %&' = 20 ()* showing the dislocation in grey, the "inefficient" PKAs in blue and the "efficient" type 1 PKAs in green and type 2 in red.

Figure 4 -

 4 Figure 4 -Unpinning probability for a shear stress = 0.95 as a function of the position of the PKA (a) along the X axis 110 ; (b) along the Y axis 1 12 ; (c) along the Z axis 11 1 . Type 1, type 2, and total unpinning probabilities are represented in green, red, and grey, respectively.

Figure 5 -

 5 Figure 5 -Unpinning probability for a shear stress = 0.975 as a function of the position of the PKA (a) along the X axis 110 ; (b) along the Y axis 1 12 ; (c) along the Z axis 11 1 . Type 1, type 2, and total unpinning probabilities are represented in green, red, and grey, respectively.

Figure 6 -

 6 Figure 6 -MD simulation results for four PKA energies at a shear stress level = 0.95 showing the dislocation in grey, the "inefficient" PKAs in blue, and the "efficient" type 1 PKA in green and type 2 in red.

Figure 7 -

 7 Figure 7 -Atoms colored by pressure and shear stress on the left and in the middle, and display of irregular atoms using the CNA tool on the right at different times after PKA creation. The corresponding videos can be found in the Supplementary data.

Figure 8 -

 8 Figure 8 -(a) Evolution of the effective volume with stress, the insert shows the dislocation using the DXA algorithm with the pinning points considered in red -(b) Pinning lifetime evolution with stress based on the distance between the pinning points in red (solid line), the hatched areas between the dashed lines represent the standard deviation, the blue line represents the experimental value obtained in Ref. [10].


Lucas Pniak

Bertrand Revol

Loïc Quéval

Jean-Sylvio Ngoua

Teu Magambo

Olivier Béthoux

Pre-sizing of a modular high power density DC/DC converter with GaN components

Keywords: Dual Active Bridge, GaN, Planar transformer, Modelling, Optimization, Aeronautical application

Finding the most appropriate architecture to achieve a function is a major engineering challenge. In the power conversion context, the question that arises is: what is the best trade-off between an architecture based on a single high-power converter and another based on the combination of multiple low-power converters? This article aims to answer this question in the specific case of the isolated interface between the low voltage (28 VDC) and high voltage (±270 VDC) buses of an aircraft. Towards this goal, the targeted power electronics converter is carefully modelled in order to obtain an accurate and fast computing model. A technological database of the different usable components is then used to feed an optimization algorithm. The execution of the latter achieves an attractive and robust result showing an excellent performance in terms of power-to-weight ratio, which is the key index of this study. The precision and speed of this computer-aided design is based on analytical models that are quick to run and therefore enable exploring almost exhaustively the search space. More specifically, the high-frequency transformer has a large relative mass and generates significant losses that are difficult to assess. A major modeling effort has been undertaken and has enable to define a simplified but convincing and accurate model (successfully assessed on a wide bandwidth using an experimental setup). It was used for the first time in this optimization process. According to the targeted aeronautical specifications, the two best solutions would make it possible to achieve a mass power density of 4.5 kW/kg, i.e. twice as high as traditional solutions. This result very clearly shows that the appropriate use of the new GaN transistors makes it possible to make a technological breakthrough. Considering the voltage and current ratings of these components, this shows that the combination of multiple partial converters is very promising and will make it possible in the future to achieve a significant increase in the compactness of electrical power conversion functions. This sizing study therefore clearly shows the potential of the standardized design of power electronics converters and the search for the best combinations (series, parallel) to meet any specific specifications. Additionnaly, the developed approach, based on a modelling effort, especially as far as the high-frequency transformer is concerned, and collection of manufacturer data, also makes it possible to limit the implementations and associated tests.

Hence, the main contributions of this article are i) the assembly of the analytical models of the various components constituting an isolated DC-DC power supply, ii) the simplification of an analytical model of the planar HF transformer, its experimental validation and its first use in a loop of optimization, iii) the fact of showing supporting figures that GaN technology enables technological breakthroughs to achieve high power density, which opens the way to modular architectures based on partial converters of electrical power, iv) the comprehensive description of a powerful pre-study design tool enabling to greatly reduce power converters' time to market.

Introduction

For reasons of energy efficiency, reliability and low cost maintenance, the aeronautical industry is gradually moving towards More Electric Aircraft (MEA), i.e. with all on-board equipment powered by electrical energy rather than pneumatic and hydraulic energies [START_REF] Tahir | More electric aircraft challenges: A study on 270 V/90 V interleaved bidirectional DC-DC converter[END_REF]. Without having fully switched to electric ancillaries, the Boeing 787 or the Airbus A380 and A350 airliners are 3 current MEA implementations. This new paradigm involves rethinking the electrical power supply architectures of the aircraft in order to ensure the best performance with respect to the key criteria of on-board mass, reliability and resilience [START_REF] Maré | Review on signal-by-wire and power-by-wire actuation for more electric aircraft[END_REF] [START_REF] Wileman | A road map for reliable power electronics for more electric aircraft[END_REF]. Specifically, the electric power architecture is based on

Preprint submitted to Mathematics and Computers in Simulation

February 24, 2023 multi on-board grids both interconnected for redundancy purpose and galvanically isolated for safety and reliability reasons [START_REF] Buticchi | On-Board Microgrids for the More Electric Aircraft-Technology Review[END_REF][START_REF] Buticchi | Multi-port DC/DC converter for the electrical power distribution system of the more electric aircraft[END_REF]. In addition, network voltage levels are quite different: 28 VDC is the standard Low Voltage DC bus (LVDC) dedicated to the control devices while 270 VDC or 540 VDC are devoted to the powerful actuators. It is also important to note that extensive considerations are aimed at further increasing the voltage of the High Voltage DC bus (HVDC): voltages of 800 V or even more are envisaged. In this context, aeronautical manufacturers have to carefully design power electronics converters [START_REF] Ounis | Multi-level integrated optimal design for power systems of more electric aircraft[END_REF], which are the key components permitting to interface LVDC and HVDC networks and accurately control the power flowing between them. This work aims to establish an effective method to accelerate the design stage.

Simultaneously with this evolution on the more electric aircraft, power electronics components diversity is significantly and rapidly evolving with the arrival at maturity of Wide Band Gap (WBG) technologies [START_REF] Trivedi | Comparative analysis of dual active bridge dc-dc converter employing Si, SiC and GaN MOSFETs for G2V and V2G operation[END_REF] such as Silicon Carbide (SiC) [START_REF] Boige | SiC power MOSFET in short-circuit operation: Electro-thermal macro-modelling combining physical and numerical approaches with circuit-type implementation[END_REF] and Gallium Nitride (GaN) [START_REF] Naradhipa | Power Density Optimization of 700 kHz GaN-Based Auxiliary Power Module for Electric Vehicles[END_REF][START_REF] Lu | A Modular-Designed Three-Phase High-Efficiency High-Power-Density EV Battery Charger Using Dual/Triple-Phase-Shift Control[END_REF]. Compared to standard silicon (Si) devices, these novel components stand higher junction temperatures enabling to decrease the heat sink size. They also faces lower losses both regarding conduction and switching ones. Such lower power losses both contributes to enhance the conversion efficiency and reduce the thermal management system size. However, the voltage and current ratings of these components do not enable them to be directly substituted for conventional ones. To achieve the expected function, it is therefore necessary to consider modular architectures based on the series and/or parallel connection of lower power elementary bricks [START_REF] Kenzelmann | Isolated DC/DC Structure Based on Modular Multilevel Converter[END_REF][START_REF] González | Modular multilevel converter for large-scale photovoltaic generation with reactive power flow and unbalanced active power extraction capabilities[END_REF][START_REF] Vidal | A multilevel modular DC-DC converter topology[END_REF], see Figure 1. This work addresses the issue of well quantifying the benefits of switching from a single converter based on directly adapted components to a multi-converter structure adopting the most appropriate components to achieve a given specification.

To address this issue, a numerical tool is developed. As technology evolves rapidly, it is vital to be able to repeat this procedure frequently in order to consider technological developments without carrying out long and costly implementations and tests. Reducing developments cost as well as time to market has indeed become a key issue. The two arising challenges are therefore to produce a simulation design tool giving relevant orders of magnitude, in order to make a relevant implementation decision, and making it possible to process a large number of parameters and possibilities in a short time.

To test the developed Simulation-Based Design tool dedicated to power electronics converters, the interconnection between LVDC and HVDC networks was chosen. This is a particularly critical case because i) the related voltage level ratio is large (around 20), ii) the link has to be power reversible and iii) voltage isolation is mandatory. The studies presented in the literature show that it is possible to reach power to weight ratios of 1 to 1.5 kW/kg by exploiting SiC semiconductors on the high voltage side and Si on the low voltage side [START_REF] Blanc | Optimization of a DC/DC dual active bridge converter for aircraft application[END_REF], about 2 kW/kg with a mixed SiC/GaN design [START_REF] Brandelero | Conception et réalisation d'un convertisseur multicellulaire DC/DC isolé pour application aéronautique[END_REF]. With a modular architecture (Figure 1), the use of only GaN semiconductors becomes possible, so that the power density can be further increased. Among the possible topologies for the elementary brick [START_REF] Forouzesh | Step-Up DC-DC Converters: A Comprehensive Review of Voltage-Boosting Techniques, Topologies, and Applications[END_REF], the Dual Active Bridge (DAB) converter, shown in Figure 2, meet the requirements of galvanic isolation, high power density, high power ratings and high efficiency [START_REF] Buticchi | On-Board Microgrids for the More Electric Aircraft-Technology Review[END_REF][START_REF] Krismer | Modeling and optimization of bidirectional dual active bridge DC-DC converter topologies[END_REF]. Several studies also show that it is suitable for series and parallel connection thanks to its current source behavior [START_REF] Teng | DCSST Multi-Modular Equalization Scheme Based on Distributed Control[END_REF][START_REF] Rolak | Efficiency Optimization of Two Dual Active Bridge Converters Operating in Parallel[END_REF]. The Dual Active Bridge was therefore selected for the elementary brick of the modular converter and comprehensively studied for the proof of concept. More precisely, its specifications are given in table 1.

When it comes to the MEA, aerospace manufacturers are facing challenges in developing smaller and lighter Rated power P to be determined between 1 and 5 kW electric power systems, to truly achieve the weight reduction benefit induced by replacing hydraulics and pneumatics with electrics. Consequently, this article adopts the power density as the only performance index. To clarify the basis of the optimization carried out, all the analytical models permitting to estimate losses and masses are detailed.

Similarly the models related to the constraints governing the power converters are specified. Additionally, it is shown that greatest modeling efforts focus on the transformer. Indeed, the high primary -secondary voltage ratio enforces the windings parallelization on the low voltage side, which greatly complicates the copper losses assessment. In this regard, an energy approach enables to provide a precise model which remains simple to implement and configure. The optimization is conducted with a Particle Swarm Optimization (PSO) algorithm [START_REF] Aubry | Optimisation du dimensionnement d'une chaîne de conversion électrique directe incluant un système de lissage de production par supercondensateurs : application au houlogénérateur SEAREV[END_REF][START_REF] Aubry | Bi-objective sizing optimization of a PM machine drive on an operating profile[END_REF][START_REF] Coello Coello | Multi-Objective Particle Swarm Optimizers: A Survey of the State-of-the-Art[END_REF]. The results presented in this article show that the GaN-based elementary brick can reach a 4.5 kW/kg power density, which is more than twice as high as what can be achieved with converters based on SiC or Si technologies (∼ 2 kW/kg). This article is organized as follows. Section 2 introduces the global scheme of the optimization process. It specifies the cost function to be optimized and the various tuning parameters. Then, section 3 describes comprehensively the power loss and mass estimation models. Section 4 presents the design constraints of the elementary brick and its components. The related models are detailed. Section 5 outlines the findings of the optimization run and shows a partial experimental validation. Section 6 draws the main conclusions and perspectives of the current work.

Simulation-Based Design

This study aims to provide a methodology for power electronics converters designers to define, by computation and before any implementation, the key parameters that their power architecture should have based solely on its specifications. The three main elements to be considered in this decisive design stage are described in this section and explained in the considered case study, namely the optimization of a high voltage ratio DAB converter. To this end, the desired purpose is first described (subsection 2.1.) to finally explain the action levers to achieve it (subsection 2.3.). Obviously, to undertake the optimization process, it is mandatory to explain the formal links that link the means to the desired end as well as the databases necessary for using these formal relationships (subsection 2.2.). Figure 3 illustrates the procedure of the method thus described.

Objective function

As aforementioned one of the remaining challenges of MEA is to achieve weight saving. This goal is a key point since mass is the main airplane fuel consumption driver while the payload is what an airline can sell. With this respect, the current study is based on the power-to-weight ratio P m as a single cost function. It is written as follows :

P m = P-ΣP i ΣM i ( 1
)
where P is the converter rated power, P i are losses in its various components, and M i are the mass of each constituent. Note that the electrical conversion efficiency is necessarily taken into account by this performance index because mass and volume of a converter are strongly determined by its losses. As a sum, the design algorithm requires on accurate models computing the losses and the masses of each elementary part.

Algorithm architecture

Solving an optimization problem requires significant computational resources, especially since the models used are numerically complex to solve and the search space is large (high number of optimization parameters and substantial variation ranges). This is indeed the case of this study. On one hand it aims providing an effective and therefore precise means of avoiding long and costly trial and error process. Secondly, it seeks to explore all the possibilities for converging towards ultra-relevant technology elementary bricks enabling the products standardization by modules combination. To achieve fast computing, the models used by the PSO optimization algorithm are analytical ones. Additionally, the considered design parameters can be either discrete or continuous variables.

The selected case study

As mentioned in the introduction section, the selected case study focuses on the most relevant DAB converter permitting to achieve a macro-converter interfacing the aircraft LVDC and HVDC buses. Its rated power and switching frequency are typical macro-parameters that have to be set by the optimization process. More specifically, eleven design parameters are considered. They could be classified in 3 categories: global operating parameters, active bridge settings and transformer characteristics. All the information are summed up in Table 2. Additionally, Figure 3 shows the general schematic of the simulation-based design tool, developed and evaluated on a DAB used for power management purposes in an embedded microgrid.

Converter modeling

The purpose of this section is to describe the collection of models needed to estimate the power to weight ratio of the Dual Active for a given set of parameters. According to the equation [START_REF] Tahir | More electric aircraft challenges: A study on 270 V/90 V interleaved bidirectional DC-DC converter[END_REF], this requires an estimate of the rated power of the DAB, the various losses it experiences and the mass of each of these components. They are detailed in the following subsections.

DAB rated power

The electrical power conversion function is performed by a combination of elementary bricks. These are based on DAB technology and controlled using a single-phase-shift (SPS) control strategy. The latter is the standard control technique consisting in driving both the LV and HV bridges by a square wave (with 50% duty-cycle) and adjusting its power p by tuning the normalized phase-shift ϕ in the range [-0.25, 0.25]. Using SPS and assuming zero losses, the power transfer from one DC bus to the other one p is determined by the following equation [START_REF] Doncker | A three-phase soft-switched high-power-density DC/DC converter for high-power applications[END_REF]:

p(ϕ) = mV p V s f L ϕ(1 -2|ϕ|) ( 2
)
where m is the transformer ratio, f is the switching frequency, V p and V s are the nominal voltage of the DC networks and L is the inductance in series with the transformer. The rated power of the Dual Active Bridge converter P is reached when p(ϕ) is maximal, corresponding to ϕ = 0.25. It gives:

P = mV p V s 8 f L . (3)
The inductor L mass and volume strongly depend both on LP 2 and also on the quantity of magnetic heat losses to be removed. The core material used as well as the switching frequency f selected highly impact these losses. Likewise, the switching frequency f deeply affects the auxiliary power derived to operate the GaN drivers. By this single equation ( 2), it could therefore already been perceived that the addressed problem of minimizing the power-toweight index is highly entangled. Subsequently, DAB modeling requires particular care. The present section and the following one deal with this demanding issue.

Loss models

3.2.1. In semiconductors Conduction losses. High Voltage (HV) and Low Voltage (LV) bridges suffer various losses. Among them, the most significant are the transistors conduction losses P c,i . They are modeled based on the on-state resistance R DS on of the GaN components, the number of components in parallel and the shape of the current i L (t) passing through the inductor L , as described in [START_REF] Turzyński | Analytical Estimation of Power Losses in a Dual Active Bridge Converter Controlled with a Single-Phase Shift Switching Scheme[END_REF][START_REF] Krismer | Modeling and optimization of bidirectional dual active bridge DC-DC converter topologies[END_REF]. The R DS on value specified at the quite high temperature of 125°C is considered in order to provide a representative insight in the toughest conditions.

P c,1 = 4 N 1 R DS on ,1 I 2 L 2 , (4)
and

P c,2 = 4 N 2 R DS on ,2 (NI L ) 2 2 , ( 5
)
where N is the transformer ratio, N 1 and N 2 the number of of GaN components in parallel per switch and I L the rms current in the inductor L.

Switching losses. The DAB converter naturally achieves Zero Voltage Switching (ZVS). It makes the switching losses P sw,i negligible, at least with regard to the HV bridge. The LV bridge is a more complicated issue. In case the parasitic inductance of the LV bridge power loop is too large, the high level of switched current (between 100 and 200 A) may induce significant losses [START_REF] Krismer | Modeling and optimization of bidirectional dual active bridge DC-DC converter topologies[END_REF]. These losses are strongly dependent on the routing technology and are therefore difficult to estimate. In the first instance, these aspects had been overlooked and the related losses neglected. It could be refined in the final design stage.

P sw,1 = P sw,2 = 0. (6)
Reverse conduction losses during dead times. During the dead time phases t dt,i , the GaN component naturally conducts the switched current I sw,i in reverse. This switched current is supposed constant as the dead time is very short compared to the switching period 1/ f . It can be calculated using the expressions provided in [START_REF] Oggier | Modulation strategy to operate the dual active bridge DC-DC converter under soft switching in the whole operating range[END_REF]. The drain-tosource voltage V ds,o f f,i drop can then be estimated using the reverse conduction characteristics of the GaN component extracted from the datasheet. It depends on the gate-to-source blocking voltage and on the switched current value. This phenomenon of reverse conduction induces important losses P dt,i which are expressed as follows. Regarding the dead time short duration, the related energy losses may be fairly low.

P dt,i = 4V ds,o f f,i I sw,i t dt,i f, i = 1, 2. (7)
Gate charge losses. At each switching period, the driver applies a voltage variation ∆V gs,i to the gate of the transistors to turn them on. It supplies a charge Q g according to the gate charge characteristic given in the datasheet. In accordance with [START_REF] Lakkas | MOSFET power losses and how they affect power-supply efficiency[END_REF], the power dissipated in this process is given by equation [START_REF] Boige | SiC power MOSFET in short-circuit operation: Electro-thermal macro-modelling combining physical and numerical approaches with circuit-type implementation[END_REF].

P g,i = 2∆V gs,i N i Q g,i f, i = 1, 2. (8)

In magnetic components

The transformer implementation considered in this study is planar technology. Thanks to their numerous advantages [START_REF] Ouyang | Overview of Planar Magnetic Technology-Fundamental Properties[END_REF][START_REF] Ngoua Teu Magambo | Planar Magnetic Components in More Electric Aircraft: Review of Technology and Key Parameters for DC-DC Power Electronic Converter[END_REF], planar transformers meet the demanding requirements of modern power electronics, especially in terms of high operating frequencies and power density [START_REF] Ouyang | Optimal Design and Tradeoff Analysis of Planar Transformer in High-Power DC-DC Converters[END_REF]. Indeed, the flat shape of its windings enables a quite homogeneous current distribution within the entire turns cross-section. More specifically this property is attractive at high frequency avoiding the impact of skin effect and proximity effects [START_REF] Mclyman | Transformer and inductor design handbook[END_REF]. Like any magnetic component, it is however impacted by two types of losses: iron losses and copper losses.

Iron losses. Iron losses are classically described using the Steinmetz formula [START_REF] Steinmetz | On the Law of Hysteresis[END_REF][START_REF] Rodriguez-Sotelo | Power Losses Models for Magnetic Cores: A Review[END_REF], P tr, f er ≃ V tr, f er k f α Bβ amp [START_REF] Naradhipa | Power Density Optimization of 700 kHz GaN-Based Auxiliary Power Module for Electric Vehicles[END_REF] with k, α and β being the Steinmetz parameters obtained from the values given in the datasheet of the magnetic materials used (least squares approximation), V tr, f er the volume of magnetic material and Bamp the maximum amplitude of the magnetic induction. Bamp is estimated using the Maxwell-Faraday Law. Considering that the LV winding composed of n 2 turns is subjected to a square wave voltage of amplitude V 2 , the induction Bamp in the transformer of core section S f is expressed as follows:

Bamp = V 2 4n 2 f S f (10)
Copper losses. When no windings parallelization is considered, Dowell's and Ferreira's models [START_REF] Dowell | Effects of eddy currents in transformer windings[END_REF][START_REF] Ferreira | Improved analytical modeling of conductive losses in magnetic components[END_REF] provide good estimates of copper losses P tr,copper . However, regarding the considered use case, the high transformer ratio causes the LV side undue hardship. To overcome the LV excessive current density and limit copper losses, parallel windings are obviously mandatory. At high frequency, the correct distribution of the currents in these windings is no longer certain and must be thoroughly assessed. To this purpose, a simplified version of the harmonic model presented in [START_REF] Pniak | AC Resistance and Leakage Inductance Estimation for Planar Transformers With Parallel-Connected Windings[END_REF] is developped in this paper. The approach is quite similar to the method proposed by Weichen in [START_REF] Chen | Model and design of PCB parallel winding for planar transformer[END_REF], but much easier to implement. It is based on a well-known physical principle called "principle of least action". Applied to electrokinetics, it states that "electric current flows along the path of least impedance" [START_REF] Feynman | The Feynman lectures on physics. Volume II : Mainly electromagnetism and matter[END_REF]. For the purpose of addressing the transformer copper issue, a short-circuit is set, for instance on the LV side of the transformer. In this case, a sinusoidal current of magnitude I 0 and pulsation ω is enforced in the primary winding, i.e. on the HV side, as represented in Figure 4. The transformer acts then as an impedance Z t composed of the windings AC resistance R ac and the windings leakage inductance L ac . It can be written as: Copper losses in transformer windings are estimated based on R ac which has to be accurately determined. The latter, as well as L ac , highly depends on the current distribution in the windings connected in parallel as it is demonstrated in [START_REF] Pniak | AC Resistance and Leakage Inductance Estimation for Planar Transformers With Parallel-Connected Windings[END_REF]. Noting I = [i 1 . . . i N ] T N×1 , the vector of the complex currents in each of the N copper layers of the transformer defined in Figure 5, R ac and L ac are the following:

Z t = R ac + jL ac ω. (11)
R ac = I † RI |I 0 | 2 , (12)
and,

L ac = I † L lk I |I 0 | 2 , (13)
Figure 5: Cross-sectional representation of a planar transformer with N conductor layers. The magnetic circuit is represented in gray, the conductors in copper color and the insulators in green.

where † is the transpose conjugate operator. According to [START_REF] Pniak | AC Resistance and Leakage Inductance Estimation for Planar Transformers With Parallel-Connected Windings[END_REF], the resistance matrix R is expressed as follows

R = L c σδh c T †                         A J -B J 2 (0) -B J 2 2A J -B J 2 . . . . . . . . . -B J 2 2A J -B J 2 (0) -B J 2 A J                         Size (N+1)×(N+1) T, (14)
with

A J (∆) = sinh (2∆) + sin (2∆) cosh (2∆) -cos (2∆) , (15)
B J (∆) = 4 cos (∆) sinh (∆) + cosh (∆) sin (∆) cosh (2∆) -cos (2∆) (16)
∆ = e c /δ. [START_REF] Krismer | Modeling and optimization of bidirectional dual active bridge DC-DC converter topologies[END_REF] and L c the mean length of one winding turns, h c the conductor layers width, σ the conductor conductivity, δ the skin depth at the pulsation ω in a material of permittivity µ 0 ,

δ = 2 µ 0 σω , (18)
and T a matrix defined as follows:

T =                  0 . . . 0 1 (0) . . . . . . 1 . . . 1                  (N+1)×(N) . ( 19
)
Referring to [START_REF] Pniak | AC Resistance and Leakage Inductance Estimation for Planar Transformers With Parallel-Connected Windings[END_REF], the leakage inductance matrix L lk is expressed as follows

L lk = T †                         L c h c σδω                         A lk -B lk 2 (0) -B lk 2 2A lk -B lk 2 . . . . . . . . . -B lk 2 2A lk -B lk 2 (0) -B lk 2 A lk                         + µ 0 L c h c                  d 1 (0) d 2 . . . (0) d N+1                                          T, (20)
where d k are insulator layers thickness defined in Figure 5, and

A lk (∆) = sinh (2∆) -sin (2∆) cosh (2∆) -cos (2∆) (21)
B lk (∆) = 4 cos (∆) sinh (∆) -cosh (∆) sin (∆) cosh (2∆) -cos (2∆) . ( 22
)
The impedance of the transformer in short circuit can then be expressed as a function of the current vector in each of the layers:

Z t = I † (R + jL lk ω)I |I 0 | 2 . ( 23
)
Note that I 0 is the magnitude of the sinusoïdal current enforced in the HV winding as described in the first paragraph of this section. According to the principle of least action, the currents flow through the path of least impedance. Consequently, the currents flowing in the parallel connected windings distribute in a way that minimizes Z t . On the contrary, the currents flowing in the turns of the high voltage winding are imposed by the source I 0 , and are constraints of the impedance minimization problem:

∀k ∈ the HV winding, i k = I 0 .

Also, assuming a perfect coupling between the two windings (no magnetizing current), the current I 0 coming into the primary winding is equal to the current going out of the secondary winding. According to the definition of the currents i k in each turn of the windings (Figure 5), this results in the sum of the N currents being zero:

N k=1 i k = 0. (25)
These M constraint linear equations on the currents (Equation ( 24) and ( 25)) can be written in matrix form defining a matrix A (size M × N) and a source current vector I s (size M × 1 containing I 0 ) :

AI = I s . ( 26
)
The explicit expressions of the matrices A and Is then depend on the distribution of the primary and secondary turns on the different layers of the planar transformer. The problem of minimizing the impedance Z t under linear constraints can then be solved using Lagrange multipliers approach. The method principle is as follows: in order to find the minimum of the function Z t (I) subjected to the constraint AI -I s = 0, the Lagrangien function L(I, λ) is defined as :

L(I, λ) = Z t (I) + AI -I s T λ. ( 27
)
The solution of the impedance minimizing problem is the stationary point of L as a function of I and λ. This means that the partial derivatives of L with respect to these variables must be zero, including the partial derivative with respect to λ. The calculation of these partial derivatives enable to express the system to be solved to determine the currents:

         R + jωL lk I + A T λ = 0 AI = I s (28)
Thus, for each set of parameters, the linear system (28) is solved to obtain the value of the currents I, which is injected into [START_REF] González | Modular multilevel converter for large-scale photovoltaic generation with reactive power flow and unbalanced active power extraction capabilities[END_REF] to estimate the ac resistance and thus the AC copper losses.

Mass estimation models

Full bridges PCBs and the related heatsink. In such pre-sizing phase, it is rather difficult to determine a full-bridge mass. Many elements are not yet well established, such as driver components, connectors, PCB surface, track size. In this context, the partial mass of these components is computed based on previous H-bridge designs. The PCB mass is assessed as a function of the number of GaN components in parallel per switch. To ensure proper heat dissipation and hence thermal management, a single heat sink is set on each of the two H-bridges. Depending on the power dissipated by the bridge and the number of transistors in parallel per switch, the minimum thermal resistance R th min between the heat-sink and the air is computed. A thorough analysis of various datasheets provided by aluminium heat-sink manufacturers permits to determine a rather invariant quantity σ m . The latter is the thermal conductance mass density.

In the case of a 4m/s forced ventilation, it equals:

σ m = 50 W/(Kkg) (29)
Subsequently, we propose the following estimation of heat-sinks masses:

M heatsink,i = 1 σ m R th min ,i , i = 1, 2. ( 30
)
Inductor. The DAB series inductor mass is also assessed based on data provided in manufacturers datasheets. A large number of high current rating inductors enables to establish another roughly invariant quantity: Φ m , the maximum flux per unit of mass. We express it as follows:

Φ m = LI sat,10% M ind ≃ 6, 7 Wb/kg ( 31
)
where I sat,10% is the saturation current corresponding to a 10% drop in the rated inductance L. Thus, an inductor, characterized by an inductance value L and a rated maximum rms current I L,max , has an estimated mass M ind defined by:

M ind = LβI L,max Φ m ( 32
)
where β is a safety factor greater than 1.0 depending on the current waveform. A β = 1.5 is considered since the maximum current carried by the inductor never exceeds 30% of its rms value I L at maximum power.

Planar transformer. The planar transformer mass M tr is related to the volume V i and mass densities ρ i of the various materials constituting the transformer. Among the components are the magnetic material, the copper and the PCB dielectric material, namely FR4. M tr is then expressed as follows :

M tr = ferrite,copper,FR4 ρ i V i . (33)
Note that the volumes V i are calculated using the geometrical parameters of the transformer.

Physical and geometric constraints

To complete these physical models and to finish the modeling of this Dual Active Bridge, it is necessary to consider the physical constraints that apply to its parameters. These constraints will define the boundaries of the optimization problem solved in this paper. Most of them concern the transformer (magnetic saturation of the core, minimal magnetizing inductance value, minimal resonant frequency or maximum winding thickness), but the thermal constraints on GaN components are also described. They are all detailed in this section.

Transformer constraints

Magnetic saturation. The induction magnitude B amp in the central core leg (having a S f section surface) must not exceed the maximum induction chosen by the program user. This maximum induction is chosen regarding the magnetic material characteristics. The constraint basic purpose is preventing any saturation in the magnetic material. In order to take into consideration the high frequency core losses, B max = 0.2 T is considered for a 3C95 material. The latter is selected for its attractive performance at the desired frequency range. Considering that a square wave of amplitude V 2 applies on the LV winding composed of n 2 turns, the inequality constraint is [START_REF] Krismer | Modeling and optimization of bidirectional dual active bridge DC-DC converter topologies[END_REF]:

B amp = V 2 4n 2 f S f ≤ B max (34)
Magnetizing inductance value. The HV magnetizing inductance L m p [START_REF] Mclyman | Transformer and inductor design handbook[END_REF] has to be larger than the DAB series inductance L by a factor k L . The latter is specified by the program user. This factor must be sufficiently large so that the reactive power exchanged with the transformer remains negligible compared to the active power exchanged through the series inductance. A choice of k L = 8 strongly limits this reactive power while enabling low magnetizing inductance values, which is necessary to fulfill the resonant frequency constraint. Thus,

L m p = A L n 2 1 1 + µ r e 0 L γ ≥ k L L, (35)
where A L is the core specific inductance (without any air gap), n 1 is the number of turns of the HV winding, e 0 is the thickness of the air gap, µ r is the core material relative permeability, and L γ is the average length of the magnetic circuit.

Resonant frequency. To guarantee the correct operation of the transformer, its open circuit resonance frequency f r must be higher than the switching frequency f of the converter. We then impose the following constraint:

f r ≥ k f f, (36)
where k f is a margin coefficient that we set equal to 3.5. To estimate the transformer resonant frequency f r , a simple equivalent model is established. As shown in Figure 6, the electrical energy stored in the windings is modeled by a simple capacitance C p brought back to the HV side. The latter mainly interacts with the magnetizing inductance L m p (large in relation to the leakage inductance L f ) around the resonant frequency which is expressed as:

f r = 1 2π L m p C p . (37)
The parasitic capacitance C p is computed from the electrical energy stored between each layer of conductors in the winding window. Three distinct cases can be distinguished, namely: two HV winding layers in interaction, two LV winding layers in interaction, and finally interaction between a HV winding layer and a LV winding layer. In the considered winding geometries, the HV winding is composed of n 1 turns while each of its conductor layers is composed of n turns On the contrary, the LV layers build a single turn; they are all connected in parallel. Noting V 1 the voltage applied to the transformer HV winding, d the thickness of the insulator of permittivity ϵ 0 ϵ r between two layers of conductor of surface S , the elementary geometrical capacitance C 0 of two turns facing each other may be expressed as:

C 0 = ϵ 0 ϵ r S nd (38)
The electrical energy stored between the conductor layers also depends on how the layers are connected to each other, as shown in Figure 7. Therefore, the energy calculation must be done for the continuous winding technique (Figure 7.a) and for the discontinuous winding technique (Figure 7.b) , as in [START_REF] Ouyang | Overview of Planar Magnetic Technology-Fundamental Properties[END_REF].

-Two HV winding layers -continuous winding. (Figure 7.a)

E pp,a = 1 2 C 0 n k=1 (2k -1) V 1 n 1 2 = 1 2 C 0 V 1 n 1 2 × n 4n 2 -1 3 ( 39
)
Figure 6: HF equivalent model of the planar transformer. R ac is the equivalent AC resistance of the HV and LV windings brought back to the HV side, L ac the leakage inductance brought back to the HV side, L m the magnetizing inductance of the HV winding and C p the parasitic capacitance brought back to the HV side representing the electrical energy stored in the transformer. -Two HV winding layers -discontinuous winding. (Figure 7.b)

Un-1 =3 V1 /n1 Un =V 1 /n1 1 0 → 2π Voltage U sup U inf = (2 n-1) V1 /n1 U1 0 → 2π 0 → 2π n-1 n U1 = ... = Un = nV 1 /n1 1 0 → 2π Voltage U sup U inf 0 → 2π 0 → 2π n-1 n U1 = VL 1 0 → 2π Voltage U sup U inf 0 → 2π 0 → 2π n-1 n U L Un-1 = VL + (n -2 )V1 /n1 Un = VL + (n -1 )V1 /n1 1 n-1 n U1 C 0 Un-1 C 0 Un C 0 1 n-1 n U1 C 0 Un-1 C 0 Un C 0 1 n-1 n U1 C 0 Un-1 C 0 Un C 0 0 → 2π 0 → 2π 0 → 2π a)
2nV 1 /n 1 V L +nV 1 /n 1 V L V 1 /n 1 0 nV 1 /n 1 nV 1 /n 1 nV 1 /n 1
E pp,b = 1 2 C 0 n k=1 n V 1 n 1 2 = 1 2 C 0 V 1 n 1 2 × n 3 (40)
-Two LV winding layers. As all the LV winding layers are electrically connected in parallel, the potential difference between turns facing each other is zero [START_REF] Ouyang | Overview of Planar Magnetic Technology-Fundamental Properties[END_REF]. It means that there is no electric field in the insulation layer and that the electrical energy stored between two LV winding layers is zero.

E ss = 0 (41)
-HV and LV winding layer. (Figure 7.c) Assuming that, for the i th HV winding layer, the left end is at potential

V L = (i -1)V 1 /n 1 , we have, E ps = 1 2 C 0 n k=1 V L + [k -1] V 1 n 1 2 = n 2 C 0        V 2 L + (n -1)V L V 1 n 1 + (n -1)(2n -1) 6 V 1 n 1 2        (42)
-Expression of the equivalent parasitic capacitance referred to the HV side.

C p = 2 E pp + E ps V 2 1 ( 43
)
PCB thickness. The PCB thickness e PCB should not exceed the winding window height e p . Hence, e PCB has to respect the following inequality constraint:

e PCB = N+1 k=1 d k + N • e c ≤ e p ( 44
)
where N is the number of conductor layers of thickness e c , and d k , k ∈ [1, N + 1], corresponds to each insulation thickness between two consecutive conductors.

Thermal dissipation. The power dissipated by the transformer must be lower than the maximum power it can naturally dissipate when it reaches the critical temperature θ c tr defined by the user. For an ambient temperature of θ amb = 20°C, a heating of 80°C is allowed, that is θ c tr = 100°C. The power dissipated by convection is assessed based on Newton's formula:

Φ c = h conv S (θ -θ amb ). ( 45
)
where h conv is the convective exchange coefficient depending both on the temperature and exchange surface orientation [START_REF]Fundamentals of heat and mass transfer[END_REF]. Finally, the heating powers that can be handled by convection is computed both regarding ferrite and transformer PCB, namely Φ c, f er and Φ c,PCB respectively. Note that at the dimensions and temperatures involved, the heat infrared radiation is negligible compared with the heat transport by convection. Thus, the following constraint must be verified,

P tr,copper + P tr, f er ≤ Φ c, f er + Φ c,PCB . (46)

H-bridges thermal constraints

According to datasheets, the GaN transistors temperature should not exceed the critical value of 150 • C. To achieve a reliable safety level, appropriate safety margins have been adopted, choosing θ j max = 120 • C as the maximum temperature threshold.

Summary

At this stage of the proceeding, all the equations describing the power electronics converter have been established based on analytical formulations. In particular, it is possible to calculate the power to weight ratio, which is the chosen objective function. As aforementioned, its expression can be written as:

F ob j = P m = P -P tr -P br,1 -P br,2 -P ind -P capa M tr + M br,1 + M br,2 + M ind + M capa ( 47
)
where P is the converter rated power, P tr is the transformer power losses, P br,1 is the HV bridge power losses, P br,2 is the LV bridge power losses, P ind is the inductor power losses, P capa is the decoupling capacitors power losses, M tr is the transformer mass, M br,1 is the HV bridge mass, M br,2 is the LV bridge mass, M ind is the inductor mass and M capa is the decoupling capacitors mass. The following section therefore describes the implementation and use of the optimization algorithm based on the models described. This part ends with the analysis of the obtained results in the specific case of a DAB converter.

Optimization implementation: results and discussion

Use Case

An algorithm exploiting these loss and mass estimation models has been developed in order to identify the most suitable converter structure for an elementary brick of a HVDC-LVDC macro-converter. As mentioned, the use case is a high ratio DAB converter. The design parameters are listed in Table 2. One of the answers that the decision-making process must provide is to determine the most suitable power for one DAB. In this study, the optimization algorithm speed benefits from the fact that some variables have limited range. For instance, galvanic isolation is achieved using a planar transformer which is a specific technology, both in terms of dedicated cores and PCB copper thickness [START_REF] Ngoua Teu Magambo | Planar Magnetic Components in More Electric Aircraft: Review of Technology and Key Parameters for DC-DC Power Electronic Converter[END_REF]. Additionally, GaN transistors are still, for the moment, under development. Consequently the number of references are quite limited. Finally, the series inductor is not custom made but selected from a manufacturer catalogue. Therefore, the overall search space is wide but could be even wider in the future with technology maturity.

Optimization problem

We look for the DAB converter sizing that minimizes the power to weight ratio P m . The optimization problem is formulated as follows: where the optimization variables X are (Table 2): the switching frequency f , the serial inductance of the DAB L, the Ferrite Core Reference FCR, the transformer winding turns thickness (i.e. thickness of the PCB conductor layer) e c , the insulation layer thicknesses in the transformer winding d pp and d ps , the air gap of the transformer core e 0 , the Winding Interleaving Configuration WIC, the HV bridge GaN reference HVGaN, the LV bridge GaN reference LVGaN, the number of transistors in parallel per switch in the HV bridge N 1 and the number of transistors in parallel per switch in the LV bridge N 2 . Note that this optimization is conducted with a specific set of constraints. The magnetic induction in the core must remain below 0.2 T. The magnetizing inductance L m p must be at least 8 times larger than the series inductance L of the DAB. The open circuit resonance frequency of the transformer f r must be at least 3.5 times greater than the DAB switching frequency f . The thickness of the transformer windings manufactured in PCB technology e PCB must be smaller than the winding window height e p . The losses in the transformer P tr must be lower than the power that can be dissipated by natural convection Φ c . And finally, the junction temperature of GaN θ j components should not exceed 120°C.

The optimization is performed by a Particle Swarm Optimization (MOPSO) algorithm [START_REF] Aubry | Optimisation du dimensionnement d'une chaîne de conversion électrique directe incluant un système de lissage de production par supercondensateurs : application au houlogénérateur SEAREV[END_REF][START_REF] Aubry | Bi-objective sizing optimization of a PM machine drive on an operating profile[END_REF][START_REF] Coello Coello | Multi-Objective Particle Swarm Optimizers: A Survey of the State-of-the-Art[END_REF] with 100 particles and 100 iterations. With a desktop computer with Intel core i7 10 th generation, 1.8 GHz and 16 GB of RAM, one single optimization takes around 10 minutes.

Results

The main results of executing the described optimization algorithm are shown in Table 3. Each of the A, B and C columns refers to the optimum sizing obtained for a given magnetic core reference. Column D describes the suboptimal achieved by the algorithm regarding E58 core reference. Obviously, one DAB sizing stands out from the others by achieving a very attractive power to weight ratio. Indeed 4.5 kW/kg is a very high level of performance twice higher than standard DAB implementation, the typical classic figure being close to 2.0 kW/kg. Besides this key piece of information, a number of further lessons can be drawn from the results analysis that bear mention.

• The best parameters set is in the middle of the search space and not at its borders, which clearly suggests that the addressed problem has a real optimal solution rather than a constrained solution.

• The best power-to-weight ratio far exceeds state-of-the-art. The optimal option (column B) achieves much higher level of performance than previous implementations do. It implies that GaN transistors are real game changers for future developments.

• The various switching frequencies are set high to limit the inductance and transformer masses. Using E64 core however induces such high LV currents as well as significant parasitic windings capacitor increase that the algorithm tends to significantly limit the switching frequency. This limitation hence negatively impacts the power-to-mass ratio. Similarly a smaller core like E43 core has a limited performance index. Indeed, lowering the LV currents mechanically reduces copper losses and enables to increase the switching frequency, which unfortunately boosts side GaN driving energy and core losses. The related global performance (column A) is significantly lower than the best sizing based on E58 core (column B).

• The sub-optimal sizing using an E58 ferrite core (column D) shows the trade-off between power density and efficiency. The increase in the number of paralleled transistors in the bridges enables to reduce the conduction losses at the cost of a higher mass. This sizing has a lower power density, but enables to gain 0.5 points of efficiency in comparison with the optimal design (column B).

• Two winding configurations stand out. The partial interleaving configurations (S-P-P-S-S-P-P-S-S-P1 (c1) and P-S-S-P-P-P-S-S-P (c2)) are the most suitable because they perform well regarding AC resistance while limiting the value of the transformer parasitic capacitance C p . This low value permits to achieve high resonant frequency and hence high switching frequency.

• Although very appropriate for high frequency operation, the planar transformer remains a key component. The optimization algorithm chooses a set of parameters maximizing the transformer RMS currents while the thermal constraint. This latter is achieved based on a global minimization of core and copper losses.

• Core and copper transformer losses are generally not balanced. Indeed, to reduce the transformer and inductor sizes, the switching frequency should be increased. In a specific frequency bandwidth, core losses may be reduced. For instance core losses of 3C95 Ferroxcube ferrite are reduced in the [100 -500]kHz bandwidth.

Conversely increasing the switching frequency may significantly boost copper losses. Skin and proximity effects indeed impact substantially these losses. However, current magnitude should be high to reach a high converter power rating. So this is why the overall optimization process favors higher copper losses than core ones.

• These very encouraging results should now be refined to consider EMC filters.

GS66516T GS66516T GS66516T GS66516T HEMT LV reference EPC2020 EPC2020 EPC2020 EPC2020 N 1 1 2 2 3 N 2 2 4 6 6

Prototype and first experimental results

The DAB candidate sizing comprehensive study offers promising prospects to achieve compact insulated DC-DC power converters. Before moving to the industrial development stage, electronic designers have to validate implementation hypotheses. In particular, dedicated research teams have to check that at least four GaN transistors may be driven using the same gate driver and routed with insignificant parasitic inductances and low resistive connections. Additionally, the carrying out should also validate the thermal management based on a single heat-sink per bridge. Regarding the planar transformer, experimental trials have to assess parasitic capacitor value and global losses at various powers.

For these specific purposes, the HV (i.e. 270V) and LV (i.e. 28V) bridges have been implemented based on GaN Systems and EPC components, respectively (Fig. 8). Planar transformers have also been developed, manufactured and tested. These elements are meant to be assembled to constitute a DAB, with a rated power of 1.9 kW and operating at a frequency of 304 kHz. Experimental basic tests have demonstrated the feasibility of i) parallelization of GaN components (4 on the LV bridge and 2 on the HV bridge), ii) bridges thermal management using a single heatsink and iii) the planar HF transformer high-quality performance, even with large current ratings. Thanks to the planar transformer prototype, we were able to validate the original model for estimating the AC resistance of the windings, presented in section 3 of this paper. Figure 9a shows the winding resistance as a function of frequency. We observe a difference of about 15% between the measurement and the model estimate at the frequency of 300kHz. This discrepancy is probably due to the 3D effects of the transformer which are not taken into account in our 2D analytical model. This is supported by Figure 9b, which shows the estimation of the winding resistance as a function of frequency by the analytical model and by 2D finite element simulation for different winding interleaving. The deviations here are very small and show that the analytical model is accurate and quick, which makes it perfectly suited for optimization.

Finally, it can be noted that the overall mass objective was successfully achieved. Table 5 compares the masses estimated by the algorithm for the optimal DAB (Table 4, column B) and the actual measured masses of the different components shown in Figure 6. Frequency (Hz) FEA simulations (asterisk markers) and the analytical model (solid lines) for the configurations P-S-P-S-P-S-P-S-P-S (fully interleaved in blue), P-S-S-P-P-S-S-P-P-S (partially interleaved 1 in yellow), and P-P-S-S-S-P-P-P-S-S (partially interleaved 2 in red).

Conclusion and perspectives

The MEA development still faces challenges to achieve aircraft auxiliary power fully based on electrical energy. Among identified issues stands the proof that switching to this form of energy is attractive in terms of reducing the mass and volume of the complete system. The present work was carried out with a view to providing an electronic power converter sizing tool enabling to accurately assess mass and losses. It is intended to be precise, rapid and generic at the same time. These specifications therefore make it possible to explore a large research domain, and therefore to carry out an exhaustive search for the best converter, i.e. the one maximizing the chosen objective function. It also makes it possible to regularly examine the optimum defined regarding the evolution of the technological characteristics of the various components used in the power converter.

For this purpose, the current study has developed analytical models of all the components implemented in the power converter, including, in particular, the HF transformer ensuring the crucial galvanic isolation function. The models could be validated during the realization of the converter obtained as an optimization process output. The case study chosen is deliberately difficult to implement and model while being particularly vital in the distribution of electrical energy within the aircraft. Based on a mono-objective optimization process (the power-to-weight ratio), the converter obtained has an attractive performance index (4.5 kW/kg, without considering EMI filters) compared to standard converters (∼ 2 kW/kg). This disruptive result is essentially due to the integral use of novel GaN components. These are connected in parallel both in the low and high voltage sides. The capability to integrate this parallelization without any extra parasitic inductance while using the same driver for each parallel switches has been fully tested on a full scale implemented prototype.

Based on these main achievements, the prospects of this work are to combine the elementary bricks so as to obtain and test the main function, i.e. LV and HV buses connection. From a methodological point of view, it seems valuable to extend the current modeling to take into account aging phenomena, making it possible to assess the function reliability as well as the required level of maintenance.

Figure 1 :

 1 Figure 1: Schematic representation of the modular converter with Input Series Output Parallel associations of the elementary brick to achieve the 540 V -28 V conversion.

Figure 2 :

 2 Figure 2: Schematic of the Dual Active Bridge converter.

Figure 3 :

 3 Figure 3: Schematic illustration of the developed method to find the most compact elementary brick enabling to achieve the desired modular macro converter.

Figure 4 :

 4 Figure 4: Equivalent circuit model of a short-circuited transformer excited by an harmonic current source I 0 .

Figure 7 :

 7 Figure 7: Left side : Possible scenarios of electrical coupling between two layers of conductors. Right side : related potential differences between the facing turns.

X

  ={ f, L, FCR, e c , d pp , d ps , e 0 , WIC, HVGaN, LVGaN, N 1 , N 2 } subject to B amp ≤ 0.2 T L m p ≥ 8L f r ≥ 3.5 f e PCB ≤ e p P tr ≤ Φ c θ j max ≤ 120 • C (48)

Figure 8 :

 8 Figure 8: (a) Transformer, (b) HV bridge, (c) LV Bridge.

Figure 9 :

 9 Figure 9: Winding resistance of the transformer as a function of frequency. (a) Comparison of the estimation of the proposed analytical model used for the DAB optimization with measurements conducted with the actual transformer presented on Figure 8. (b) Comparison between 2-DFEA simulations (asterisk markers) and the analytical model (solid lines) for the configurations P-S-P-S-P-S-P-S-P-S (fully interleaved in blue), P-S-S-P-P-S-S-P-P-S (partially interleaved 1 in yellow), and P-P-S-S-S-P-P-P-S-S (partially interleaved 2 in red).

Table 1 :

 1 Dual Active Bridge converter specifications

	Parameter	Nominal value Range of variation
	Input voltage V e	270 V	[210; 330] V
	Output voltage V s	28 V	[16; 32] V

Table 2 :

 2 Design parameters

	Category	Symbol Parameter
	DAB operating parameters	f L	Switching frequency DAB inductance value
		FCR	Ferrite Core Reference
		e c	Conductor thickness
	Geometric parameters	d pp	insulator thickness between conductor of the same winding
	of the transformer	d ps	insulator thickness between conductor of different winding
		e 0	air gap thickness
		WIC	Winding Interleaving Configuration
		HVGaN HV bridge GaN reference
	Parameters related	LVGaN	LV bridge GaN reference
	to the bridges	N 1	Number of transistors in parallel in the HV bridge
		N 2	Number of transistors in parallel in the LV bridge

Table 3 :

 3 Parameter sets maximizing the power mass density of the elementary conversion brick.

	Parameters		Optimal Solutions		Sub-optimal solution
	P m (kW/kg)	3.78	4.47	4.1	4.08
	P (kW)	1.01	1.92	2.4	1.99
	η (%)	97.46	97.7	97.9	98.2
	f (kHz)	410	304	176	272
	L (µH)	22.8	16.2	22.4	17.5
	core reference	E43/10/28-3C95 E58/11/38-3C95 E64/10/50-3C95	E58/11/38-3C95
	Winding configuration	(c2)	(c1)	(c1)	(c1)
	Resonant frequency (kHz)	1448	1090	632	1020
	P tr,core (W)	4.4	3.64	2.2	3.9
	P tr,copper (W)	3	8.0	11.3	7.8
	e (µm)	140	175	210	210
	d ps (µm)	470	535	275	470
	d pp = d ss (µm)	325	230	200	250
	e 0 (µm)	125	210	270	190
	HEMT HV reference

Table 4 :

 4 Mass estimation and comparison.

	Component	Estimated	Measured
	LV Bridge	115 g	130 g
	HV Bridge	96 g	100 g
	Transformer	185 g	180 g
	Inductance	37 g	50 g
	TOTAL	433 g	460 g

This notation means that the winding (c1) is composed of 10 layers of conductor. Layers 1, 4, 5, 8 and 9 are for the secondary winding, and layers

2, 3, 6, 7 and 10 are for the primary winding. The numbers correspond to those defined in Figure5.