domain: phys chim spi

text:
Christopher Glaser
email: christopher.glaser@onera.fr

Jouke Hijlkema

J Ér Ôme Anthoine

APPROXIMATING IDEALIZED HYBRID ROCKET FUEL PORT GEOMETRIES USING STEPS

Keywords: Hybrid propulsion, stepped geometry, profile optimization, regression rate

ABSTRACT:

In this work, a novel approach to increase the performance of a Hybrid Rocket Engine (HRE) by approximating any idealized fuel port geometry by steps inside the fuel grain is proposed. This work serves as a first investigation into the effect of single steps inside the fuel grains of a typical HRE. In numerical simulations, it is shown that forward-and backwardfacing steps have the potential to increase the average regression rate as they promote the turbulence and heat transfer inside the turbulent boundary layer. Both types of steps are more effective at lower oxidizer mass fluxes. Finally, as a first estimation, a set of best practices for the distribution of steps along a profile could be derived.

INTRODUCTION

The unique features of Hybrid Rocket Engines (HREs) such as safety and simplicity make them promising candidates for low cost space access. The downsides like low regression rates, slower transients, possibly high fuel slivers and low maturity, however, need to be addressed.

In short, we propose to approximate an idealized fuel port profile by a set of cylinders with different inner diameters (refer to Fig. 1). The flow will therefore display flow phenomena resembling a flow over backward-(BFS) and forward-facing steps (FFS). This way, the engine can profit both from the regression rate enhancement due to the steps, and exhibits the benefits of an idealized profile. Adding to this, the design is easy to manufacture, since there is no need for additive manufacturing of complex port geometries.

In order to elaborate the idea, this article is structured as follows: First, the existing research regarding regression rate enhancement in HREs through steps and diaphragms is briefly collected and analysed. Diaphragms are considered here, because the flow phenomena they induce lie somewhere between a forward-and backward-facing step. The literature review is followed by a numerical study using the multi-physics solver CEDRE of ONERA. Here, the effect of a single step on the flow inside a hybrid rocket combustion chamber is investigated. Furthermore, in evaluating several step configurations such as forward-and backward-facing steps with different heights, an important height-to-length relationship for the recirculation zone can be derived. An exten-sive case study for a HRE with steps (and its unique interdependence of fuel addition from the fuel surface, mixing, heat transfer and regression rate) is lacking to the knowledge of the authors.

THE CONCEPT

Starting point is the fact that the characteristics of HREs can be altered by changing the fuel port profiles. These profiles can be optimized according to different criteria. For example, a genetic algorithm (previously developed at ONERA) [START_REF] Hijlkema | A presentation of a complete design cycle for optimised hybrid rocket motors[END_REF] can be fed with different optimization criteria (e.g. constant oxidizer-to-fuel ratio or low fuel residuals) to calculate an idealized fuel port profile. However, producing these geometries with conventional fuel casting methods is difficult. Hence, additive manufacturing is often the only possibility to create the desired profiles. Moving to larger scales, this leads to increased complexity and cost. This is where a second observation comes into play. Over the years, several researchers investigated the effect of diaphragms and other obstacles (such as steps) in the HRE combustion chamber. Overall, the results were the same; due to the turbulence and enhanced mixing after the obstacles (introduced by the recirculation zone) the heat transfer and the regression rate increased.

The proposed design approach is therefore the following: If an idealized profile can be approximated by a set of steps (as depicted in Fig. 1), one could theoretically achieve three major advantages: 1) The profile no longer needs to be printed but can be approximated by a finite number of fuel cylinders with different inner diameters. 2) The idealized profile provides the characteristics it was optimized for, e.g. lower O/F shift during operation. 3) Due to the turbulence introduced by the steps, the regression rate and mixing increases. In this paper, the main goal is to derive a set of considerations and guidances on how to best distribute a set of steps along any desired profile. This is done by simulating the impact of single BFS and FFS for different step heights on the fuel surface temperature, regression rate, fuel mass flow and recirculation zone length (RZL). These data allow to understand which is the most efficient way to distribute a set of BFS and FFS along a predefined profile.

DIAPHRAGMS AND STEPS IN THE LITERATURE

The concept of steps or diaphragms in the fuel port to promote turbulence and mixing is not new. ON-ERA employed a diaphragm in the LEX (Lithergol Experimental) rocket in the 1960s [START_REF] Duban | La fus ée sonde lex[END_REF], making it the first diaphragm in a hybrid rocket motor with flight heritage. Gany and Timnat [START_REF] Gany | Parametric study of a hybrid rocket motor[END_REF] observed average regression rate increases of 50 % when employing a diaphragm in the middle of their polyester fuel grain already in 1972. Grosse [START_REF] Grosse | Effect of a diaphragm on performance and fuel regression of a laboratory scale hybrid rocket motor using nitrous oxide and paraffin[END_REF] extensively conducted studies on diaphragm form and position in paraffin HREs (a work that has been continued at University of Padua [START_REF] Bellomo | Investigation of effect of diaphragms on the efficiency of hybrid rockets[END_REF]). The common results show regression rate increases downstream the diaphragms by up to 90 % for lab-scale motors and combustion efficiency increases from 76 % to 95 %. Zhang et al. [START_REF] Zhang | Numerical simulations on unsteady operation processes of N2O/HTPB hybrid rocket motor with/without diaphragm[END_REF] pointed out that the regression rate after the diaphragm peaks consistently at the point of the reattachment of the flow. This effect was exploited by Kumar and Kumar [START_REF] Kumar | Effect of diaphragms on regression rate in hybrid rocket motors[END_REF] (numerically) and Dinesh and Kumar [START_REF] Dinesh | Utility of Multiprotrusion as the Performance Enhancer in Hybrid Rocket Motor[END_REF] (experimentally): Both teams placed multiple diaphragms in the fuel grain to further extend the effect of diaphragms. They evaluated that the diaphragms are best to be placed after the recirculation zone of the preceding diaphragm. Quantitatively speaking, for Kumar and Kumar the ideal spacing was 8-10 times the height of the diaphragm in numerical simulations [START_REF] Kumar | Effect of diaphragms on regression rate in hybrid rocket motors[END_REF], which corresponds to the observed RZL. Interestingly, Dinesh and Kumar observed the optimal spacing to be about 2.7 times of the RZL. Concerning the use of steps in HREs, first records of a hybrid rocket motor with a backward facing step can be found in Korting et al. [START_REF] Korting | Advanced hybrid rocket motor experiments[END_REF] in 1987. It was shown that a BFS at the port entrance (hence, technically the step is not inside the fuel grain) doubled the local regression rate after the step. Noteworthy, the peak of the regression rate enhancement appeared at the point of the re-attachment of the flow after the step, as has been pointed out for diaphragms also. Lee et al. [START_REF] Lee | Combustion instability for hybrid rocket motors with a diaphragm[END_REF] proposed the use of a backward facing step instead of a diaphragm to avoid the pressure oscillations often induced by diaphragms. The stepped geometry was achieved by combining two fuel grains with different inner di-ameters. With a 7.5 mm step, the regression rate increased by about 50 % downstream the backward facing step. Sakashi et al. [START_REF] Sakashi | Effectiveness of Concave-convex Surface Grain for Hybrid Rocket Combustion[END_REF] used multiple steps with alternating inner diameters. Their so-called concave-convex design effectively resembles a sequence of BFS and FFS with a maximal height difference of 9 mm. In doing so, the regression rate could be increased experimentally by up to 100 %. Kumar and Joshi [START_REF] Kumar | Regression rate study of cylindrical stepped fuel grain of hybrid rocket[END_REF] used 5 cylindrical fuel grains with either 25 mm or 35 mm inner diameter. Combining them, a BFS/FFS fuel port could be created that increased regression rates between 30-50 %. Contrary to the multi-diaphragm research, for steps, the optimal spacing due to the RZL has not been exploited.

Distributing a set of BFS and FFS along an optimized profile has neither been proposed nor investigated. This is where our proposal to distribute a set of BFS and FFS to follow a predefined optimized profile comes into play.

NUMERICAL SET-UP

The simulations are carried out with the multiphysics solver CEDRE of ONERA [START_REF] Refloch | Cedre software[END_REF]. This section provides a brief overview of the underlying numerical models and assumptions of the hybrid rocket simulations. As this paper is focused on the application of the model rather than the numerical implementation, the concepts are only briefly described. For a detailed and thorough understanding it is highly recommended to refer to the work of Durand et al. [START_REF] Durand | Turbulence Modeling Effects on Fuel Regression Rate in Hybrid Rocket Numerical Simulations[END_REF][START_REF] Durand | D éveloppement et validation exp érimentale d'une mod élisation num érique pour la simulation d'un moteur-fus ée hybride[END_REF] who implemented the HRE regression rate model into CEDRE. The base-line of the simulation is based on the HYCAT (Hybride avec CATalyseur ) facilities of ONERA. The HYCAT family uses H 2 O 2 and High-Density Polyethylene (HDPE) as propellants. The hydrogen peroxide is decomposed with a catalyst into hot water steam and gaseous oxygen. This is why only the gas phase has to be considered for the simulations, decreasing also the numerical complexity. As for the pyrolysis products of HDPE, only gaseous ethylene (C 2 H 4 ) is assumed. Together with the decomposed hydrogen peroxide this leads to the following chemical species in the simulations:

H 2 O, CO, CO 2 , O 2 and C 2 H 4 .
Given the significant difference of fuel regression rate (mm/s) and oxidizer flow velocity (m/s), the flowfield is considered as quasi-steady. The simulations are 2D-axisymmetric and the flow is described by the compressible Reynolds-Averaged Navier-Stokes equations. The set of equations is closed with the Shear Stress Transport (SST) model of Menter [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF], combining the benefits of the k-ϵ model in the core flow and the k-ω model near the walls. Due to the fuel injection at the fuel surface, the correction of Wilcox [START_REF] Wilcox | Reassessment of the scaledetermining equation for advanced turbulence models[END_REF] for the specific dissipation at the fuel wall is considered. The fuel inlet is modeled with an Gas Surface Interaction (GSI) model using an Arrhenius-type description to calcu-late the regression rate based on the fuel surface temperature.

Chemical reaction

The combustion inside an HRE is described by a turbulent diffusion flame. The two-step kinetic chemical reaction denotes according to Westbrook and Dryer to [START_REF] Westbrook | Simplified reaction mechanisms for the oxidation of hydrocarbon fuels in flames[END_REF]:

C 2 H 4 + 2O 2 -→ 2CO + 2H 2 O (1) CO + 0.5O 2 + H 2 O ←→ CO 2 + H 2 O (2)
The detailed rates of production and destruction can be found in Ref. [START_REF] Durand | Turbulence Modeling Effects on Fuel Regression Rate in Hybrid Rocket Numerical Simulations[END_REF].

Gas-surface interaction model

The GSI model at the fuel surface used in CEDRE of Durand et al. [START_REF] Durand | Turbulence Modeling Effects on Fuel Regression Rate in Hybrid Rocket Numerical Simulations[END_REF][START_REF] Durand | D éveloppement et validation exp érimentale d'une mod élisation num érique pour la simulation d'un moteur-fus ée hybride[END_REF] is based on an energy-and mass balance at the fuel surface. It allows to couple the fuel surface temperature to the regression rate, allowing to calculate the regression rate rather than fixing it at the fuel surface.

According to mass conservation at the fuel surface, the mass flux of the fuel (G fuel ) can be written as:

G fuel = ρ s • ṙ , (3)
with ρ s being the solid fuel density and ṙ the regression rate. The production rates ωi of each species i are considered source terms and yield the mass balance:

y i • G fuel -ρD m,i ∂y i ∂η = ωi , ( 4
)
with D m being the diffusion coefficient, y the mass fraction and η the radial coordinate. We consider the sum of the mass source terms ωi to be equal to the total fuel mass flux and ethylene as the sole pyrolysis product:

N i=1 ωi = G fuel . ( 5
)
Moreover, the energy balance (gas phase left, solid phase right) denotes to :

λ g ∂T ∂η + Q rad - N i=1 ωi • h i = λ s ∂T ∂η s -G fuel • h fuel . ( 6
)
λ g/s : Thermal conductivity gas/solid T : Temperature Q rad : Net radiation h : Enthalpy N : Total number of species Assuming a quasi-steady heat flux through the solid fuel:

λ s ∂T ∂η s = G fuel • c p,fuel • ∆T , ( 7
)
Eq. 6 simplifies to (by neglecting radiative heat transfer) [START_REF] Durand | Turbulence Modeling Effects on Fuel Regression Rate in Hybrid Rocket Numerical Simulations[END_REF]:

λ g ∂T ∂η = N i=1 ωi • h i -G fuel • ∆H • f,fuel , (8)
with c p,fuel and ∆H • f,fuel being the heat capacity and formation enthalpy respectively. Finally the production rate of ethylene (C 2 H 4 ) -hence the fuel mass flow G fuel -is given by the Arrhenius law:

G fuel = A pyr • exp (- E a R • T ) . ( 9
)
The values for the pre-exponent A pyr , the activation energy E a and the formation enthalpy ∆H • f,fuel for HDPE are listed in Tab. 1.

Table 1: Properties of HDPE [14]

A pyr [ kg m 2 s ] E a [ J mol ] ∆H • f,fuel [ J kg ] 3.5*10 6
125.6*10 5 -6.365*10 5 ± 1.14*10 4

Computational domain

The computational domain in this article is based on the HYCAT facilities at ONERA, allowing to be able to validate the simulations and their continuation in the future. In order to decrease the number of influences on the results of the fuel surface, at this stage of the research the HYCAT family is approximated by a 2D axisymmetrical academic domain, as depicted in Fig. 2 for the BFS case. The open source software GMSH is used to create the meshes. The total fuel grain length is 240 mm. The inlet fixes the oxidizer mass flux (G ox ) and the inlet temperature. The fuel surface is modeled according to the GSI presented in Sec. 4.2 and the outlet is a pressure outlet which fixes the pressure to a value representative to a nozzle. The mesh is chosen to be Cartesian because it lessens the number of cells (as compared to triangles) and allows to refine the mesh size close to the fuel surface to approximate the boundary layer region. The minimal cell size at the fuel surface is calculated to ensure an y+ < 1 for all possible total cell counts. The inner diameter before the step is fixed to 25 mm as this represents the minimum diameter of HYCAT. The diameter after the step is dependent on the step height (h). For a maximum step size of 10 mm this leads to a diameter after the step of 45 mm. Changing from BFS to FFS domain is conveniently achieved by switching the inlet and outlet conditions. This signifies that for the FFS cases, the diameter after the step is always fixed to 25 mm and the diameter before the step changes according to the step height. Translated to a real world application this means that a cylindrical fuel grain with two different inner diameters is simply turned around.

Mesh convergence

Since at this stage of the research we are not yet able to validate the simulations experimentally, it is of upmost importance to prove that the results are independent from the mesh. In order to do so, six different meshes with increasing cell count (from roughly 15,000 to over 400,000) are created for a BFS domain with a step height of 2 mm. All meshes ensure a y+ < 1 at the fuel surface. The results are presented in Tab. 2. Mesh 1 (14,902 cells) serves as starting point and the cell count is approximately doubled with each mesh. After each individual simulation is converged, the average values of surface pressure, surface temperature, total fuel mass flow and RZL are gathered for 10,000 iterations with 100 data points for the same physical time. Mesh 1 presents the base case, the values for meshes 2-6 give the relative change compared to the preceding simulation. For example, for mesh 4, the total fuel mass flow increased by 0.16 % as compared to the fuel mass flow of mesh 3. The boundary conditions of the mesh convergence study are loosely based on the conditions inside a typical HYCAT firing and can be found in Tab. 3.

Looking at Tab. 2, it becomes obvious that the average values of pressure and temperature converge very fast and could be considered mesh independent starting from mesh 1. For the fuel mass flow, however, the values can be considered convergent somewhere between mesh 2 and 3. The largest changes in values can be seen for the RZL.

After the first refinement to mesh 2, the length of the recirculation zone increased by almost 20 %. This is due to the fact that the coarse mesh 1 is not able to properly capture the eddies that are separating after the step. After mesh 3, the RZL finds a more stable value. However, in order to correctly judge the mesh independence for the RZL, it is necessary to also consider the shape of the recirculation zone. Fig. 3 depicts the appearance of the recirculation zone for Meshes 1,4 and 5. It becomes apparent that the RZL for case 1 is not at all developed and too short. For mesh 4, however, the length and the height of the RZL seems to have a more consistent appearance. Nonetheless, a smaller second recirculation zone inside the bigger zone is observable, which could have a non-negligible effect on the simulations. Moreover, the contour of the flame zone is less consistent as compared to mesh 5. Mesh 6 only changes the RZL by under 2 %, but uses double the cells. For this reason, mesh 5 with 208,307 cells is considered as the most suitable mesh for the numerical test campaign in terms of mesh independence and computational cost. Before going into the details of the results, we present in short the basic characteristics of the flow phenomena of a flow over a BFS and FFS respectively. Fig. 4 shows the temperature fields and velocity streamlines for a representative flow without step, a BFS with 5 mm height and a FFS with 5 mm height. For the BFS case, the flame zone detaches after the step and re-attaches after the recirculation zone. The flow before the step is unaltered. For the FFS case, however, the step already lifts up the flame zone the closer it gets to the step and a first recirculation zone is already formed before the step. After the step, the flow re-attaches considerably faster than for the BFS case. This is due the fact that in the BFS case, the diameter increases and therefore the flow velocity decreases for a constant mass flow. In the case of the FFS, on the other hand, the diameter decreases and the flow accelerates. This fundamental difference has to be kept in mind when discussing the flow phenomena in the succeeding sections.

Test matrix

The simulations will be carried out on the domain displayed in Fig. 2. The boundary conditions are being held close to the HYCAT data. Using only a single step at this stage allows to single out the effects of the steps. The interactions between multiple steps will be considered in future work. The following major investigations need to be carried out:

Height study:

The impact on different step heights h.

Position study:

The effect of the axial position of the step (x/l).

Mass flow study:

The importance of the oxidizer mass flux G ox = ṁ A given by the mass flow ṁ. These studies are carried out both for the BFS and FFS cases. Between the three different parametric campaigns, all other values are being held constant. The different inlet conditions, step heights and positions can be found in Tab. 5, the boundary conditions (based on typical HYCAT firings) that are similar throughout all calculations are listed in Tab. 4. The results of all studies are collected and their impact on the average regression rate, recirculation zone and temperature and regression rate profiles assessed in the respective subsections.

Backward facing step cases

In the following section, first the BFS cases are analyzed separately.

Average regression rate

Fig. 5 depicts the spatially and temporal averaged regression rate values for the BFS cases limited to behind the step. It allows to evaluate the impact of all parametric studies in one single image. The green graph represents the reference cases without any step. The blue line and the crosses correspond to the average regression rate values downstream the step. The label S2 i.e. corresponds to a step height of 2 mm, the label S5 end refers to the case of 5 mm step height at the position x/l = 0.75. Additionally, the Marxman [START_REF] Marxman | Turbulent boundary layer combustion in the hybrid rocket[END_REF] fit is given as

ṙ = a • G n ox ( 10
)
where a and n are parameters used to best fit the data points. For the Marxman fit, only the data points of the mass flow study are considered. This allows to extrapolate approximated values even outside of ranges that have been simulated.

Comparing the Marxman fits of the reference cases and the BFS cases, the most prominent observation is that the average regression rate of the BFS cases is inferior at mass fluxes above approximately 100 kg/m 2 s. Below this value, the average Figure 5: Regression rate after step for BFS regression rate behind the step is higher than the reference cases without step. Moreover, the cases S1, S2, and S5 end are below the reference cases; whereas the cases S10 and S5 front are above the reference cases. This observation hints at the existence of two competing effects: Directly after the BFS, the flame zone is further from the fuel surface, therefore decreasing the fuel surface temperature. On the other hand, the step introduces a recirculation zone that promotes mixing and heat transfer. If the recirculation zone is too small (as it seems to be the case for 1 mm and 2 mm for higher mass fluxes), the regression rate after the steps never recovers from the flame zone being further from the fuel surface directly after the step. This is also why the case of S5 front, where the step is positioned at the beginning of the fuel port rather than the middle, is superior to the references cases: The turbulence induced by the step has more time and space to compensate the decrease of the surface temperature directly after the step. In order to discuss this effect we introduce the temperature and regression rate surface profiles of the different cases in the next part.

Surface profiles

Fig. 6 shows the temperature distribution along the fuel port surface along the x-axis of the fuel port. Moreover, the gray boxes indicate the length of the recirculation zone induced by the step. The green dotted line specifies the profile of the respective reference cases. Keeping in mind that the diameter of the fuel port increases after the step, we need to consider a different reference case after the step than before. For example, with a mass flow of 100 g/s in a 5 mm case, the reference G ox before the step yields 203 kg/m 2 s. After the step it decreases to around 104 kg/m 2 s due to the larger diameter. The step position is at 0 m distance.

With the help of these surface profiles, the aforementioned competing effects of recirculation zone and flame distance from the surface becomes clearer: In all cases, the surface temperature directly after the step drops significantly, given that the flame is further from the surface. However, the surface temperature increases exactly until the reattachment point (after the recirculation zone). From there on, the heat transfer and the mixing is increased, clearly visible in steeper temperature profiles as compared to the reference case. With the help of Fig. 6 it can be explained why the cases for S1 and S2 are inferior to the reference cases: The increased heat transfer induced by the step is too small to compensate the loss of surface temperature directly after the step. Starting with S5, the augmented heat transfer is considerably higher than the loss of temperature immediately after the step and can in total increase the surface temperature. Fig. 7 displays the regression rate profiles of Figure 6: Temperature of the BFS height study the same cases as Fig. 6. The regression rate profiles follow the exact same trend as the temperature profiles. This is little surprising given that the surface temperature was modeled in the GSI (recall also Sec. 4.2) to be the main influence of the regression rate. Knowing this, we will continue from here on to only consult the regression rate profiles for the sake of efficiency. The regression rate profiles confirm the competing effects: Directly after the step, the regression rate descends well below the reference value. At the re-attachment point, the regression rate peaks (an effect also observed by Refs. [START_REF] Gany | Parametric study of a hybrid rocket motor[END_REF][START_REF] Zhang | Numerical simulations on unsteady operation processes of N2O/HTPB hybrid rocket motor with/without diaphragm[END_REF][START_REF] Korting | Advanced hybrid rocket motor experiments[END_REF]) and after a smaller rebound continues with a steeper profile than the reference cases.

The effect of the steeper profile can be further illus-Figure 7: Regression rate of the BFS height study trated with the results of the step position study in Fig. 8. Here, the step height is constant at 5 mm but the position moved from 25 % (front) to 75 % (end) of the total grain length. It becomes apparent that a BFS at 75 % of the grain length diminishes the effect of the step. Contrarily, for a BFS at 25 %, the steeper regression rate profile considerably increases the regression rate for over the half of the total fuel grain. This explains also why the spatially averaged regression rate of the S5 front case in Fig. 5 is superior to the other S5 cases.

Figure 8: Regression rate profile for different step positions. The step height is constant at 5 mm.

Regression rate for averaged mass flux

To conclude the investigation of the BFS simulations, we draw another comparison to the reference cases using the spatial and temporal averaged regression rate. This time, however, we also average the mass flux for the BFS cases. As known, the mass flux G ox is given by the mass flow ṁ and area A as:

G ox = ṁ A = ṁ 0.25 • π • D 2 . ( 11
)
For stepped cases like the BFS and FFS cases, the diameter changes after the step. That is why we define an average mass flux G ox,avg as:

G ox,avg = ṁ A avg , ( 12
)
where the average area A avg depending on the step position (x/l) is defined as:

A avg = 0.25 • π • [(x/l) • D 2 1 + (1 -(x/l)) • D 2 2 ] (13)
with D 1 and D 2 being the diameter before and after the step respectively. With this expression we can now compare the temporal and spatial averaged regression rate of the whole domain with the references case. Fig. 9 illustrates the results. The overall trend is similar to 5, nonetheless, the average mass flux where the BFS cases become superior to the reference moves up to around 130 kg/m 2 s. It becomes apparent that lower mass fluxes and step heights of 5 mm and higher yield the best results, while the position of the step should be as early as possible.

Figure 9: Regression rate of averaged mass flux

Forward facing step cases

Following the same procedure as for the BFS cases, this section discusses the outcome of the FFS simulations.

Average regression rate

Fig. 10 reports the average regression rate behind the FFS in red and compares them to the reference cases without a step. The prominent observation is that -contrary to the BFS cases -all FFS cases exhibit a higher regression rate than the reference cases. However, the difference declines with increasing mass flux. Preliminary, this leads to the conclusion that FFS are better suited to increase the regression rate through turbulence and enhanced heat transfer than BFS. This conclusion will be elaborated further in the next section when looking at the surface profiles. The reason why most of the FFS cases are clustered at the same G ox value is that for the FFS domain, the diameter after the step is fixed to 25 mm, whereas the diameter before the step changes depending on the step height. That is why in Fig. 11 a more detailed view of the area labeled Detail A is provided.

Here, it is to be noted that the S1 case has about the same regression rate as the reference cases. Starting with S2, the difference grows. Interestingly, there seems to be a negligible difference between the S5 and cases. A trend opposite to the BFS cases where the BFS was by far superior to all other cases. Moreover, for the BFS cases, the regression rate was highest when the step position is in the front. For the FFS, on the other hand, the step position at the end exhibits the highest regression rate. This observation will also be further investigated with the help of the surface profiles in the next section.

Surface profiles

Before taking the surface profiles into account, it is reminded that we will concentrate solely on the regression rate profiles. This decision is justified, since we showed in the BFS section that the temperature and the regression rate profiles follow the exact same trend, as it is expected from the implementation of the Arrhenius law in the GSI. Likewise to the BFS surface profiles, Fig. 12 compares the regression rate to the reference cases while indicating also the RZLs. Due to the diameter being smaller after the step (leading to increased G ox ) the reference cases before and after the step are different in order to refer to the proper mass flux.

The first observation is that contrary to the BFS cases, the influence of the FFS is already noticeable before the step. This is due to the step 'pushing' the flame further away from the surface the closer it gets to the step. As soon as the recirculation zone before the step is reached, the regression rate drops dras-tically. This leads to an overall smaller regression rate before the step as compared to the reference cases. For the 10 mm case, this effect is the most visible. After the step, the regression rates show Figure 12: Regression rate of the FFS step height study a steep increase with a distinct peak exactly at the re-attachment point. For the BFS cases, this peak was visible, but less pronounced. Moreover, for the BFS cases, the recirculation zone is considerably larger than for the FFS cases. This translates to the regression rate peak being reached earlier than for the BFS cases. After the peak, the regression rate follows approximately the same angle as the reference case. Interestingly, the difference between the 5 mm and 10 mm step is almost negligible after the step, whereas the negative effects before the step itself are considerably higher. This leads to the conclusion that the the FFS should not be higher than 5 mm. Another fundamental differences to the findings of the BFS cases, where the 10 mm case is to be preferred.

If we turn now to the step position study that is detailed in Fig. 13, we can observe another discrepancy to the BFS cases: For the FFS cases, the step position further downstream increases the regression rate the most. At the 25 % position, the peak is not only smaller but also less pronounced. Additionally, the decrease in regression rate before the step is noticeable already before the recirculation zone starts. Given that the turbulent boundary layer and its flame zone is less developed at the 25 % position, we can assume that the disturbance of the flow at the earlier position is stronger.

Figure 13: Regression rate profile for different step positions. The step height is constant at 5 mm.

Regression rate for averaged mass flux

Following Eq. 12 that was derived to average the mass flux not only after the step but over the whole domain, Fig. 14 displays the averaged regression rate over the average mass flux. In line with Fig. 10, nearly all FFS cases are superior to the reference cases. Additionally, the effect increases with decreasing mass flux. Unexpectedly, the S5 end case drops below the expected trend of the FFS Marxman fit, although in Fig. 13 this case seemed superior. Since the step in this case is at 75 % of the total canal length, the increased regression rate after the FFS only contributes to 25 % of the total canal length. Therefore the disturbances at 75 % of the total length caused before the step dominates during the averaging process. Nonetheless, the S5 end case lies still above the reference cases.

Height-to-length relationship

In this section, we put the length of the recirculation zone (x) in relation to the height of the step (h) that is inducing the recirculation zone. This yields the height-to-length relationship x/h. In order to compare the different relations amongst each other, we introduce the step Reynolds number Re h , following a common approach in the literature [START_REF] Chen | A review of Backward-Facing Step (BFS) flow mechanisms, heat transfer and control[END_REF][START_REF] Awasthi | Pressure fluctuations produced by forward steps immersed in a turbulent boundary layer[END_REF]:

Re h = U • ρ • h µ , ( 14
)
Figure 14: Regression rate of averaged mass flux with U being the free stream velocity, ρ the flow density and µ the dynamic viscosity. All values are calculated at the domain inlet. Fig. 15 displays the RZL for the BFS case. For the FFS, the two recirculation zones (before and after the step) are given. Additionally, a power law fit is plotted for each data-set and the coefficient of determination (R 2 ) specified. The values of the BFS recirculation zones lie between 4.9 and 6.9 and their evolution is coherent with the trends that are found in (non-HRE) literature of BFS flows (eg. Ref. [START_REF] Chen | A review of Backward-Facing Step (BFS) flow mechanisms, heat transfer and control[END_REF]): With increasing Re h , the x/h ratio increases and shows saturating behaviour. This also can serve as one possible explanation why the BFS cases benefit from lower mass flows and thus, lower Re h values. The regression rate for BFS cases has shown to be inferior to the reference cases within the recirculation zone. It is typically after the recirculation zones that the regression rate overtakes the reference cases (refer to Fig. 7). Hence, with smaller recirculation zone relative to the step height, the BFS cases become more competitive. Given that the x/h decreases with decreasing mass flow, this could be an explanation for the trend. For the x/h ratio before the FFS, the trend also is well in line with what is reported in (non-HRE) FFS literature [START_REF] Awasthi | Pressure fluctuations produced by forward steps immersed in a turbulent boundary layer[END_REF]: With increasing mass flow and Reynolds number, the recirculation zone ratio before the FFS step decreases. Unexpectedly, the data for the x/h ratio after the FFS is inverse to what is found in (non-HRE) literature [START_REF] Awasthi | Pressure fluctuations produced by forward steps immersed in a turbulent boundary layer[END_REF]. Literature reports the x/h after the FFS to increase with the Reynolds number. However, for our FFS cases, the ratio decreases with increasing Reynolds number. The reason for this behaviour could not be fully understood yet and is part of on-going investigations in the future. One possible explanation can be the differences between the cases in the literature and our HRE related set-up. For example, in the FFS literature, usually the flow is not axisymmetric but on flat surfaces. In the HRE FFS case, we have an axisymmetric geometry, where the flow significantly accelerates after the FFS due to the smaller diam-Figure 15: Length-to-height ratio of recirculation zones eter. Moreover, according to the Marxman law, the regression rate increases with increasing mass flux after the FFS. This leads to increased blowing from the surface, especially at the end of the recirculation zone where a pronounced peak of regression rate was observable. All these alterations to the typically planar literature case could serve as explanation for the discrepancy of the (non-HRE) literature and our results.

CONCLUSION AND OUTLOOK

In this work, we proposed a novel approach to increase the performance of HREs. The concept is to approximate any idealized profile inside a fuel port with a set of BFS and FFS (recall also Fig. 1). This way, potentially, the engine can profit from the idealized fuel profile, while also exhibiting higher regression rates due to the turbulence inducing steps. Since there is few available data on the effects of BFS and FFS in HREs, the first step to promote the idea was to understand the effect of single steps on the regression rate. These results will build the foundation of future algorithms to distribute multiple steps. Hence, the following main observations are to be recalled:

1. FFS cases exhibit in all but one case higher regression rates than the non-step references.

2. BFS cases are assumed to have a 'breakeven' point (below around 130 kg/m 2 s oxidizer mass flux) where they are superior to the reference cases.

3. Both BFS and FFS have shown to be more effective with decreasing oxidizer mass flux.

4. In all cases, the peak in regression rate was observed to be at the re-attachment point.

5. For BFS, the higher the step, the better the regression rate enhancement; for FFS, the augmentation has shown to be saturated for 5 mm step height.

6. For BFS cases, the regression rate increases more after the recirculation zone than within, that is why a step at the beginning of the fuel port yielded better results than at the other positions downstream. This trend is considerably less pronounced for FFS, and might even be inverse.

Concluding these investigations, a first careful conclusion can be drawn, serving as an initial best practice for the distribution of steps along an approximated profile: It is advised to use fewer but taller and longer BFS, preferably upstream the fuel port, while the height of the FFS should be limited to 5 mm. FFS are less susceptible to its length and axial position inside the fuel port.

For the future, it is necessary to validate these first findings not only numerically but also experimentally. Moreover, this work has been concentrated solely on single steps. It is of importance to confirm the results for multiple steps and validate if the interaction of several steps is benefiting the regression rate increase or maybe even diminishing it. Another concern is the evolution of the steps with time. If the steps dissolve rapidly as the burn continues, the turbulence enhancing effects would also get lost. Nonetheless, the profile that was once approximated by steps would be now smooth and could still benefit from the profile optimization (while it was less complex to manufacture as compared to printed solutions). This will be investigated in continuation of this work.

In conclusion, we have demonstrated that BFS and FFS have a turbulence and regression rate enhancing effect on HREs and can thus be potentially used to approximate fuel port geometries. The results of the parametric study will serve as a starting point to develop an algorithm to distribute steps along an idealized profile. Even though the research is still in an early stage, we believe that BFS and FFS can be a promising addition to HREs, especially when they are used to minimize manufacturing effort while also increasing the performance.

ρ

  . . . . . . . . . . . . Density ω . . . . . . . . . . . . Production rate η . . . . . . . . . . . . Radial coordinate λ . . . . . . . . . . . . Thermal conductivity µ . . . . . . . . . . . . Dynamic viscosity A . . . . . . . . . . . . Area Apyr . . . . . . . . . . Arrhenius pre-exponent cp . . . . . . . . . . . Heat capacity D . . . . . . . . . . . . Diameter Dm . . . . . . . . . . . Diffusion coefficient Ea . . . . . . . . . . . Activation energy G . . . . . . . . . . . . Mass flux ∆H • f . . . . . . . . . . Formation enthalpy h . . . . . . . . . . . . Step height hi . . . . . . . . . . . Enthalpy I turb . . . . . . . . . . Turbulent intensity N . . . . . . . . . . . Total number of species P . . . . . . . . . . . . Pressure Q rad . . . . . . . . . . Radiation ṙ . . . . . . . . . . . . Regression rate R 2 . . . . . . . . . . . Coefficient of determination Re . . . . . . . . . . . Reynolds Number R . . . . . . . . . . . . Gas constant T . . . . . . . . . . . . Temperature U . . . . . . . . . . . . Velocity x/h . . . . . . . . . . . Length-to-height ratio x/l . . . . . . . . . . . Step position y . . . . . . . . . . . . Mass fraction

Figure 1 :

 1 Figure 1: The proposed stepped concept

Figure 2 :

 2 Figure 2: Computational domain for BFS case

Figure 3 :

 3 Figure 3: Recirculation zone after the step for different meshes

Figure 4 :

 4 Figure 4: Temperature contours of axisymmetric HRE simulation

Figure 10 :Figure 11 :

 1011 Figure 10: Regression rate after step for FFS

Table 2 :

 2 Mesh convergence study Mesh no. Cell count Pressure [Pa] Temperature [K] Fuel mass flow [g/s] Recirculation zone [mm]

	1	14,902	3,602,121	924.74	5.978	8.56
	2	26,199	+0.025 %	+0.139 %	+2.11 %	+19.77 %
	3	56,263	+0.027 %	+0.032 %	+0.56 %	+2.91 %
	4	96,683	-0.0002 %	+0.011 %	+0.16 %	+3.77 %
	5	208,307	+0.013 %	+0.051 %	+0.9 %	-2.73 %
	6	412,097	-0.004 %	+0.017 %	+0.6 %	+1.87 %

Table 3 :

 3 Boundary conditions of mesh convergence

	T inlet	P outlet	G ox	y H2O y O2 I turb
	1000 K 36 bar 200 kg m 2 s	0.59 0.41 1 %

Table 4 :

 4 Constant boundary conditions

	T inlet	P outlet	y H2O,inlet y O2,inlet I turb
	900 K 36 bar 0.59	0.41	11 %

Table 5 :

 5 Values for parametric studies

	Study	h [mm]	x/l [-]	ṁ [g/s]	Labels
	Height	[1, 2, 5, 10]	0.5	100	[S1, S2, S10]
	Position	5	[0.25, 0.5, 0.75]	100	[S5 front, S5 end]
	Mass	5	0.5	[50 , 100, 200]	none

ACKNOWLEDGEMENTS

The authors would like to thank Dr. Levard for his guidance on CEDRE and the meshing of the computational domain.

The project leading to this application has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 860956. It is part of the ASCenSIon project, an Innovative Training Network (ITN) to advance space access capabilities.


Armand Barbot

Riccardo Gatti

Unsupervised learning for structure detection in plastically deformed crystals

Keywords: Structure Detection, Unsupervised learning, Molecular Dynamics, Crystal, Plastic deformation

Detecting structures at the particle scale within plastically deformed crystalline materials allows a better understanding of the occurring phenomena. While previous approaches mostly relied on applying hand-chosen criteria on different local parameters, these approaches could only detect already known structures.

We introduce an unsupervised learning algorithm to automatically detect structures within a crystal under plastic deformation. This approach is based on a study developed for structural detection on colloidal materials. This algorithm has the advantage of being computationally fast and easy to implement. We show that by using local parameters based on bond-angle distributions, we are able to detect more structures and with a higher degree of precision than traditional hand-made criteria.

Introduction

Molecular dynamics simulation is a powerful method allowing to simulate at the particle scale different materials such as colloidal systems [START_REF] Boattini | Unsupervised learning for local structure detection in colloidal systems[END_REF], glassy materials [START_REF] Barbot | Rejuvenation and shear banding in model amorphous solids[END_REF][START_REF] Barbot | Characterization of shear bands and plasticity in model glasses at the atomic scale[END_REF] or metallic nanocrystals [START_REF] Amodeo | Mechanical properties and dislocation nucleation in nanocrystals with blunt edges[END_REF]. To help in interpreting the simulations results, being able to determine the local structure at the particle-scale is essential. To do so, several approaches were developed, mainly relying on local order parameters to describe the surrounding environment of each particle (number of neighbours, angles formed with the neighbours, ...) to detect underlying substructures in the simulated atomistic sample. Among these methods, we can cite the Bond Orientational Order (BOO) parameter [START_REF] Steinhardt | Bond-orientational order in liquids and glasses[END_REF], the Common Neighbours Analysis (CNA) [START_REF] Honeycutt | Molecular dynamics study of melting and freezing of small Lennard-Jones clusters[END_REF], or the Bond-Angle Distribution (BAD) [START_REF] Ackland | Applications of local crystal structure measures in experiment and simulation[END_REF]. Such methods were applied with success to study several phenomena such as crystal nucleation [START_REF] Sanz | Out-of-equilibrium processes in suspensions of oppositely charged colloids: liquid-to-crystal nucleation and gel formation[END_REF], melting [START_REF] Noori | Study of the effect of grain size on melting temperature of Al nanocrystals by molecular dynamics simulation[END_REF] or plasticity [START_REF] Ackland | Applications of local crystal structure measures in experiment and simulation[END_REF][START_REF] Stukowski | Structure identification methods for atomistic simulations of crystalline materials[END_REF][START_REF] Amodeo | Atomistic Simulations of Compression Tests on Ni 3 Al Nanocubes[END_REF]. However, they are mostly relying on hand chosen criteria and thus only works for already known structures.

Recently, different approaches using machine learning to detect substructures from local order parameters were developed. The first models, which relied on supervised learning [START_REF] Geiger | Neural networks for local structure detection in polymorphic systems[END_REF][START_REF] Dietz | Machine-learning approach for local classification of crystalline structures in multiphase systems[END_REF][START_REF] Boattini | Neuralnetwork-based order parameters for classification of binary hard-sphere crystal structures[END_REF], suffered from the same problem as the hand-chosen techniques: they could only be trained to find an expected structure in the studied systems. However, a recent work from Boattini et al [START_REF] Boattini | Unsupervised learning for local structure detection in colloidal systems[END_REF] successfully applied an unsupervised learning method using the BOO parameter to automatically detect different structures within a colloidal material without relying on a priori knowledge of the underlying structures. This method also had the advantage of being easy to implement and computationally fast.

While the BOO parameter is able to discriminate between the different structures in crystalline materials, it does not perform well when the simulated system is under deformation [START_REF] Lechner | Accurate determination of crystal structures based on averaged local bond order parameters[END_REF]. Being able to automatically detect substructures within crystal under plastic deformation, including some previously undetected and unexpected case, would open the way to a better understanding of crystal plasticity at the atomic scale. This requires to find a local order parameter suitable to discriminate substructures at the atomistic level for elastically and plastically deformed systems.

In this study, we present a method inspired by the paper from Boattini et al [START_REF] Boattini | Unsupervised learning for local structure detection in colloidal systems[END_REF] to automatically study and detect the different substructures within a crystal under plastic deformation appearing at the atomistic scale. This approach relies on the BAD parameters designed for single crystals to describe the environment around each atom, as highlighted in [START_REF] Ackland | Applications of local crystal structure measures in experiment and simulation[END_REF]. This local parameter was already used in previous structure detection approaches relying on bond-angle distributions, showing the ability to discriminate structures associated with crystal plasticity [START_REF] Ackland | Applications of local crystal structure measures in experiment and simulation[END_REF][START_REF] Amodeo | Atomistic Simulations of Compression Tests on Ni 3 Al Nanocubes[END_REF].

Our method consists on extracting the most pertinent BAD parameters using an autoencoder neural networks [START_REF] Rumelhart | Learning internal representations by error propagation[END_REF][START_REF] Dietz | Machine-learning approach for local classification of crystalline structures in multiphase systems[END_REF] before applying two clustering models: K-means [START_REF] Lloyd | Least squares quantization in PCM[END_REF]Mac-Queen, 1967) and DBSCAN [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF]. It has the advantage of being computationally fast and easy to implement. We finally test our approach on a FCC single crystal of Nickel under plastic deformation.

Methods

This section is divided in three main parts. We first describe the system on which we apply our structure detection method. We then depict the BAD parameters we use to capture the local environment of each particle. The last section focuses on the algorithms we use to detect the local atomistic substructures in the system from the BAD parameters.

System

In this study, we perform Molecular Dynamics (MD) simulations using LAMMPS [START_REF] Thompson | LAMMPS -a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales[END_REF] to plastically deform a Ni FCC defect free single crystal. These simulations are then used as benchmark to apply our method, based on unsupervised learning, developed to detect the different structures emerging in the system during plastic deformation.

We first create a cubic sample with edge lengths equal to 20nm and free surfaces oriented along the 100 directions as shown on Fig 1 (a). The interatomic interaction linking the atoms is the EAM potential for Ni [START_REF] Adams | Self-diffusion and impurity diffusion of fcc metals using the five-frequency model and the Embedded Atom Method[END_REF]. The system is then equilibrated for 5 ps at 5K using Nose-Hoover thermostat [START_REF] Nose | A unified formulation of the constant temperature molecular dynamics methods[END_REF][START_REF] Hoover | Canonical dynamics: Equilibrium phasespace distributions[END_REF]. A deformation is imposed using uniaxial compression with a flat indenter along the [100] direction, with a strain rate of 10 8 s -1 . The flat indenter is modeled as an infinite plane exerting a repulsive force on atoms defined by

F (r) = -K(r -R) 2 , ( 1
)
with r being the atom position, R the plane position and K the force constant.

Here we choose K=1000 eV/ Å3 . The temperature is maintained at 5K during the compression through the Nose-Hoover thermostat.

To determine when the first plastic event occurs during the deformation, we first measure on the fly the compressive stress applied on the system by the indenter using a method similar to the one used in [START_REF] Amodeo | Atomistic Simulations of Compression Tests on Ni 3 Al Nanocubes[END_REF] and looking for the first stress drop corresponding to the first plastic event. The stress-strain curve obtained from the compression is shown on Fig 1 (b), with the first plastic event represented by a vertical red line.

In our simulations, we observe the nucleation of Shockley head partial dislocations with burger vectors:

{ 1 6 [-1 1 2], 1 6 [1 -1 2], 1 6 [-1 -1 2], 1 6 [-1 - 1 -2]}.
The dislocation lines are thus accompanied by a stacking fault in the HCP structure [START_REF] Hull | Chapter 5 -Dislocations in Face-centered Cubic Metals[END_REF].

Bond Angle Distribution Parameter

To describe the local structural environment around each particle, we measure the Bond-Angle Distribution (BAD) parameters for each of them, based on the method described in [START_REF] Ackland | Applications of local crystal structure measures in experiment and simulation[END_REF]. For each atom i, we first extract its nearest neighbors N b (i) using the adaptive cutoff method as described in [START_REF] Stukowski | Structure identification methods for atomistic simulations of crystalline materials[END_REF]. This method has the advantage of being parameter-free, so can be directly applied to any structure, while being just slightly slower than using a fixed cutoff [START_REF] Stukowski | Structure identification methods for atomistic simulations of crystalline materials[END_REF].

We define the bond-angle θ jik with j, k ∈ {1..N b (i)} being two nearest neighbors of i. We then estimate {cos(θ jik )} for all the N b (i)(N b (i) -1)/2 neighbor pairs of atom i. From it, we describe nine ranges of {cos(θ jik )} over which the BAD parameters {χ l } l∈{0..8} are estimated as shown in Table 1. These ranges are based on those from [START_REF] Ackland | Applications of local crystal structure measures in experiment and simulation[END_REF] which were optimised to differentiate crystal structures. The idea here is to count the number of bond-angle cosines in each ranges to obtain the value of the corresponding {χ l } l∈{0..8} .

Table 1: Definition of the nine bond angle cosines ranges over which the BAD parameters χ are estimated. cos(θ jik ) is the cosines of the angle between r ij and r ik , with j, k being two nearest neighbors of the atom i. The ranges were taken from [START_REF] Ackland | Applications of local crystal structure measures in experiment and simulation[END_REF] For instance, let us consider a case where a atom i has six bond-angle cosines. If two of them are in the range [-1.0, -0.945] and four of them are in the range [-0.945, -0.915], then for this atom, we will have χ 0 = 2, χ 1 = 4 and {χ l } l∈{2..8} = 0.

Note that in this paper, the BAD parameters were calculated using the package Pyscal [START_REF] Menon | pyscal: A python module for structural analysis of atomic environments[END_REF].

Unsupervised learning

From the knowledge of the BAD parameters for each atom, the aim of our method is to determine automatically the different structures present in the system. This part is done through unsupervised learning and can be decomposed on two main tasks. First we use an autoencoder to determine the good number of parameters to describe the local structures in our system. Then, we apply a perturbative method to the autoencoder to determine which parameters are the most effective. From this, we focus on the most effective parameters applying clustering and classification methods to extract the different substructures of our system.

Note that as the surface atoms can be easily filtered from their number of nearest neighbours as in [START_REF] Amodeo | Atomistic Simulations of Compression Tests on Ni 3 Al Nanocubes[END_REF]. , i.e. atoms for which N b < 11, we filter out the surface atoms using this criteria and focus on the structures within the inner atoms.

Autoencoder

As explained above, the aim of this section is to extract the number of pertinent parameters to determine the different structures of our system. Indeed, out of the nine BAD parameters, some could be irrelevant or redundant to describe the local structural environment around each atom. This task is here performed using an autoencoder based neural network to perform dimensionality reduction [START_REF] Goodfellow | Deep Learning[END_REF][START_REF] Chen | Collective variable discovery and enhanced sampling using autoencoders: Innovations in network architecture and error function design[END_REF][START_REF] Boattini | Unsupervised learning for local structure detection in colloidal systems[END_REF].

The autoencoder network can be divided in two parts: the encoder part and the decoder part as shown on Fig 2 . The encoder aims at encoding the input into a lower dimension: the bottleneck. Then the decoder will try to reconstruct the input from the bottleneck. The output of the network will then be the reconstruction of the input done by the decoder from the bottleneck. In practice, the input of the encoder is

χ(i) = {χ l } l∈{0..8} (i) ∈ N 9 , with i ∈ [1..N ],
N being the number of atoms over which the training in performed.

We train the encoder to reduce the difference between the input and the output for a bottleneck dimensions going from 1 to 9.

To determine the pertinent dimension of input, we measure for each bottleneck N bneck dimension the mean square error between the input χ and the output χ as

E N bneck = 1 N N i=1 χ(i) -χ(i) . (2)
We then look from which dimension we have the minimum difference between the input and the output.

For more details about the implementation of the autoencoder, please refer to the paper [START_REF] Boattini | Unsupervised learning for local structure detection in colloidal systems[END_REF] which served as a reference for this section.

Relevant parameters identification

After determining the relevant number of parameters to study the substructures in our system, we want to go beyond the black box of the neural network and determine which parameters are most pertinent for our study. To do so, we use the input perturbation method [START_REF] Yao | Forecasting and analysis of marketing data using neural networks[END_REF][START_REF] Scardi | Developing an empirical model of phytoplankton primary production: a neural network case study[END_REF][START_REF] Gevrey | Review and comparison of methods to study the contribution of variables in artificial neural network models[END_REF][START_REF] Olden | An accurate comparison of methods for quantifying variable importance in artificial neural networks using simulated data[END_REF].

The principle of this method is to take the trained autoencoder neural network with the optimal bottleneck dimension, and to perturbate one of the input parameters by increasing its value with a fixed amount. In this work, we increase the chosen input value by 10% which was shown to give satisfactory results [START_REF] Boattini | Unsupervised learning for local structure detection in colloidal systems[END_REF].

After perturbing one of the inputs, we look at how much this perturbation influenced the output by calculating:

E pert k = 1 N N i=1 χ(i) -χpert k (i) , (3)
with k ∈ [0..8] being the index of the perturbed input, χpert k the output after the perturbation and χ the input without perturbation.

If the perturbed input k was considered during the training as carrying no pertinent information to reconstruct the input at the output, the perturbation will not propagate to the output. This will lead to a very small value of E pert k . On the contrary, if the perturbed input k is necessary to obtain an output similar to the input, its perturbation will lead to a higher value of

E pert k .
From this, we determine the relative importance RI k of the input k as:

RI k = ∆E k 8 j=0 ∆E j , (4)
with ∆E k = E pert k -E being the difference between the autoencoder mean square error after perturbation of the input k and without perturbation. From this relative importance index, we can determine the importance of each input during the autoencoder training.

Clustering and classification

In the last section, we described how to isolate the most relevant BAD parameters to study the local structure. In our system most of the atoms will be in two structures, FCC and HCP (in the stacking fault), which corresponds to close values of BAD parameters [START_REF] Ackland | Applications of local crystal structure measures in experiment and simulation[END_REF]. While most structure detection approaches focus on the already known main crystalline structures, i.e. FCC, BCC and HCP [START_REF] Stukowski | Structure identification methods for atomistic simulations of crystalline materials[END_REF], the goal of our method is also to be able to automatically detect substructures without having an a priori knowledge about them.

To this end, we chose to apply two clustering methods on the relevant BAD parameters space, one to detect the main structures containing most of the atoms, and the second one to focus on the substructures containing few atoms but located far from each other in the parameter space. The clusters extracted by these two clustering approaches will then be combined to identify the different substructures of the studied system.

In the next sections, first we introduce the K-means clustering method [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF][START_REF] Lloyd | Least squares quantization in PCM[END_REF] we use to detect the main structures. Then, we detail the Density-Based Spatial Clustering of Applications with Noise (DBSCAN) [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF] which allows to detect efficiently more isolated clusters.

While K-means is able to efficiently determine the cluster of a new data point without having to re-apply the clustering over the whole data set, it is not the case for DBSCAN. We thus need to apply a classification algorithm over the clusters detected by DBSCAN. We thus present in the last part the logistic regression classifier we use to perform this task.

Kmeans

In order to detect the main substructures present in the system, we use the Kmeans clustering algorithm [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF][START_REF] Lloyd | Least squares quantization in PCM[END_REF] as implemented in sci-kit learn [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF]. The principle of the method is the following: let us consider a set of N data points in the parameter space that we want to separate into K clusters. We name {x i } i∈[1..N ] the position of the data point in the parameter space and {µ j } j∈[1..K] the center of each cluster, commonly called "centroid".

The aim of the K-means algorithm is to choose optimal centroids in order to minimize the inertia defined by

IN (K) = N i=1 min j∈[1..K] ( x i -µ j 2 ). ( 5
)
If the points are located in clusters, the inertia will be minimized by placing the centroids at the center of each cluster.

While the K-means requires the number of clusters as an input, the optimal number of cluster can be obtained with the following procedure: we first apply K-means to our data set for varying number of clusters K and calculate the corresponding inertia IN (K). In this paper, we did it for K ∈ [1..10]. Then we apply the elbow method to the function IN (K) following the protocol described in [START_REF] Salvador | Determining the number of clusters/segments in hierarchical clustering/segmentation algorithms[END_REF]. From this method, we can obtain the value of K from which the inertia stop decreasing with K, thus corresponding to the optimal number of clusters for K-means clustering.

After training the K-means algorithm over the data set for the optimal cluster number K, we obtain the centroids locations {µ j } j∈[1..K] of the different clusters. From this information, we can determine in which cluster is a new data point by looking at its closest centroid. This allows to classify new data depending on their cluster without having to retrain the algorithm each time.

The K-means method being only looking at the distance between the data points and the centroid, it does not perform well for elongated clusters orfor clusters with irregular shapes. However, it performs very well if most of the data points are located in the same position in space. As the main crystalline structures are each associated with a specific position in the BAD parameter space, the K-means algorithm is well suited to detect the clusters associated with these structures.

DB Scan

To be able to detect clusters with smaller number of points and non-isotripic shape, we apply another clustering method: the Density-Based Spatial Clustering of Applications with Noise (DBSCAN) [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF] as implemented in sci-kit learn [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF].

The principle of the algorithm is the following: we input two parameters: (1) the maximum distance between two points l max and (2) the minimum number of points per cluster N min . In this method, two data points located at a distance l < l max are considered as neighbors. From it, the algorithm will class the data points into three categories:

(i) if a data point has at least N min -1 neighbors, it will then be considered as a "core point".

(ii) if a data point has less than N min -1 neighbors but at least one of them is a core point, it will be considered as a "border point".

(iii) the remaining data points will be considered as "noise".

The core points and the border points are then regrouped together to form clusters, located at a distance l cluster > l max from each other. In our system, many atoms can be associated with the same values of BAD parameters. As the DBSCAN only look at the distance between points, having many points located at the same spatial position will only slow down the clustering without improving the result. We thus apply the DBSCAN over the list of the unique spatial positions of our data set. Finally, the BAD parameters being integers, we use for the distance parameter the value l max = 1. We also determined that N min = 5 gives satisfactory results for our data set.

As this method is not influenced by the number of point at a given spatial position but only by the distance between the unique spatial positions in the data set, it would be able detect clusters associated with structure containing fewer atoms but isolated from the clusters corresponding to the main crystalline structures.

Logistic regression classifier

In the opposite of K-means, if we want to know from DBSCAN in which cluster is a new data point, we would need to re-perform the clustering over the whole data set, plus the new data point. But this solution would be very inefficient and time consuming. One solution to solve this problem is to train a classification algorithm over the clusters detected by DBSCAN.

The classification algorithm takes as an input the training data set and as the output the cluster found by DBSCAN. It will then decompose the parameter space into region; each associated with a cluster. When a new data point is added, it will be attributed a cluster depending on its location in the parameter space.

In this study, we use the logistic regression classifier [START_REF] Brzezinski | Logistic regression modeling for contextbased classification[END_REF] implemented in scikit-learn as it simple, fast and gives satisfactory result when applied to the clusters which detected by DBSCAN in this system. More details about this classifier can be found in sci-kit learn [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF].

Cluster combination

In the previous parts, we described two clustering methods. One to detect the main structures in the system: the K-means algorithm, the other one to detect the substructures containing few points, more isolated in the parameter space: the DBSCAN. Then, both K-means and the association of DBSCAN with the logistic regression classifier will separate the parameter space into regions corresponding to detected clusters. Thus, each data point will be associated with two clusters: the one detected by K-means, the other detected by DBSCAN.

In order to have a precise detection of the structures present in the system, we combine the results obtained from the two clustering methods. Finally, we obtain N Kmeans • N DBSCAN clusters, with N Kmeans and N DBSCAN corresponding to the number of clusters found by K-means and DBSCAN, respectively.

Results and discussions

In this section, we present the results obtained by following the procedure described in the previous sections. First we analyse the data obtained we get from the autoencoder: the relevant number of the BAD order parameter necessary to detect the substructures in our system as well as which of the χ(i) = {χ l } l∈{0..8} (i) ∈ N 9 parameters are the most significant. We then apply the clustering methods on the extracted parameter, later analysing their performances. Finally, we examine the different substructures detected by this method, comparing the outcomes with other structure detection methods.

Autoencoder output analysis

After compressing our Ni crystal until the first plastic event occurs, as described in Section 2.1, we record at a regular time step interval the atomic positions of the simulated nano-object. We then measure BAD parameters for each atoms as detailed in Section 2.2, filtering out the surface atoms as explained in Section 2.3. It turned out that computing the BAD parameters of a single time step is sufficient to train the autoencoder and to observe the appearance of different clusters. Here, the time step used to train our algorithm corresponds to an applied strain of ∆ = 0.05% after the first plastic event.

Fig 3 shows the mean square error between the autoencoder input and output for different bottleneck dimensions. We can observe that the error reaches a plateau for a bottleneck dimension of N = 3. We thus consider here that the pertinent number of BAD parameters to study the structure in our system is therefore three. Our next step is to determine which are the three most pertinent BAD parameters out of nine we computed for each atom. To this end, we calculate the relative importance (RI) of each parameter applying the perturbation method described in Section 2.3.2. The results are shown in Fig 4 . From this figure, we can clearly observe the most significant three parameters: the triplet {χ 4 , χ 5 , χ 7 }. From here on in this study, we will only focus on these three BAD parameters to study the substructures within the studied Ni nano-object.

Clustering methods combination

We determined that the most pertinent BAD parameters to study the emerging substructures in our numerical sample are {χ 4 , χ 5 , χ 7 }. Here, we apply two clustering algorithms on these parameters, K-means and DB-SCAN, and then we combine the outcome of the obtained clusters. The dataset used for the clustering processes is the same that we used for extract BAD parameters using the autoencoder (we also verified that using different atomistic configurations extracted after the onset of plasticity gives similar results) .

On Fig 5 (a), we show a scatter plot representing for each atom their position in the {χ 4 , χ 5 , χ 7 } space, as well as the result of the K-means clustering. Note that as the BAD parameters are integers, one position in the {χ 4 , χ 5 , χ 7 } space can correspond to many atoms. Each point represented here corresponds to at least one atom. On this figure, we can see that the algorithm detected two clusters: one in dark blue and the other in yellow.

As detailed in Section 2.3.3.1, the aim of using the K-means clustering is to detect the main crystalline structures, located close to each other in the BAD parameter space [START_REF] Ackland | Applications of local crystal structure measures in experiment and simulation[END_REF] From this figure, we can see that the large majority of the atoms are concentrated on two close positions in the {χ 4 , χ 5 , χ 7 } space: {21, 13, 24} and {24, 12, 24}. On [START_REF] Ackland | Applications of local crystal structure measures in experiment and simulation[END_REF], these two configurations were defined as corresponding to HCP and FCC, respectively. As HCP corresponds to the stacking fault structure in a FCC crystal [START_REF] Hull | Chapter 5 -Dislocations in Face-centered Cubic Metals[END_REF], we can expect that most of the atoms will be in these two structures. From what we can see, K-means seems to have only succeeded in separating these two configurations located within the central block. This will be confirmed in the next part focusing on analysing the structures detected by the clustering algorithms.

While the K-means algorithm successfully separated the two locations in the {χ 4 , χ 5 , χ 7 } space, we can see on From this figure, we can see that the DBSCAN successfully managed to separate into three clusters the three blocks observed in the scatter plot. After measuring the proportion of data points contained in each cluster, it appears that the central cluster (in purple) contains 99.8% of the points, while the two remaining clusters contains 0.1% each. This confirms that most of the atoms are located in the main structures, corresponding to the central block, making it difficult for the K-means clustering to detect the two smaller blocks.

Finally, to detect the main crystalline structures while being able to discern other structures containing fewer atoms in the system, we combine the results obtained by the two clustering methods (K-means shown on

Analysis of detected structures

On the previous section, by combining two clustering methods, we managed to extract six different clusters from our dataset by a projection in the {χ 4 , χ 5 , χ 7 } space. We now focus on analyzing the structures associated with each cluster. To this end, we show on Fig 7 (a) a snapshot made with OVITO [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with OVITO -the Open Visualization Tool[END_REF] of the Ni nano-cube deformed plastically, with an applied strain of ∆ = 0.02% beyond the yield point. On this snapshot, the atoms are colored corresponding to their clusters. For instance, the atoms contained in the yellow cluster on From this figure, the atoms colored in green and red follows what looks like a dislocation line, while the atoms in dark blue seems to correspond to the stacking fault created after the glide of a partial dislocation. By using the algorithm DXA [START_REF] Stukowski | Automated identification and indexing of dislocations in crystal interfaces[END_REF] implemented in OVITO [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with OVITO -the Open Visualization Tool[END_REF], we could confirm our prediction: the green and the red atoms are ). On panel (a), we observe that most of the atoms belongs to the yellow cluster which we associate with the FCC structure. On (b), we show the same snapshot after filtering out the FCC structure. We can observe that the combination of the two clustering methods is capable to separate the atoms linked with dislocation lines (green and red atoms) and those correlated with stacking faults (dark blue atoms).

surrounding a partial dislocation line. The relative position of the green and red atoms depends on the burger vector as well as on the direction of the dislocation line (see the appendix Appendix A for more details). From this, we can deduce that the dark blue atoms correspond to HCP stacking fault structures generated by the gliding of a partial dislocations (as mentioned above in FCC crystals stacking faults have a HCP structure).

The two remaining clusters detected in in Fig 6 c), the orange and the sky blue one, are very rarely observed during the time steps recorded during the deformation. The orange structure is observed when two dislocation lines are located close to each other while the sky blue one is detected during the interaction between a dislocation line and a stacking fault, as shown later.

On Fig 8, we compare our structure detection method, shown in panel (a), with three other existing methods: (i) hand-chosen criterion using the BAD parameters from [START_REF] Ackland | Applications of local crystal structure measures in experiment and simulation[END_REF] (Fig 8 (b)), (ii) the CNA method [START_REF] Honeycutt | Molecular dynamics study of melting and freezing of small Lennard-Jones clusters[END_REF][START_REF] Stukowski | Structure identification methods for atomistic simulations of crystalline materials[END_REF] implemented in OVITO ,which is one of the most commonly used method for structure detection (Fig 8 (c)) and which consists of using hand-chosen criterion on a radial-distribution based local order parameter, and (iii) the autoencoder approach using the BOO parameters from [START_REF] Boattini | Unsupervised learning for local structure detection in colloidal systems[END_REF] which served as an inspiration for this paper (Fig 8 (d)). For this comparison, we filtered beforehand the surface atoms following the method explained in section 2.3.

On [START_REF] Ackland | Applications of local crystal structure measures in experiment and simulation[END_REF], (c) the CNA method [START_REF] Honeycutt | Molecular dynamics study of melting and freezing of small Lennard-Jones clusters[END_REF] from OVITO [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with OVITO -the Open Visualization Tool[END_REF], and (d) the combination of BOO and autoencoder as detailed in [START_REF] Boattini | Unsupervised learning for local structure detection in colloidal systems[END_REF]. Note that for the four structure detection methods shown here, we removed the surface atoms beforehand.

Out of the three methods used for comparison, two were mostly designed to detect the main crystalline structures (FCC, HCP or BCC): the handchosen criterion using the BAD parameters and CNA. We can still remark that CNA, which also relies on hand-chosen criterion, is able to capture the location of the dislocation line but with less accuracy than our method: the number of atoms close to the line defect is much larger than in our method. Also, this method is not able to capture sub-structures associated with plasticity. The atoms outside of the main crystalline structures, such as those in the line defect, are labelled as "other" (in red in Fig 8 (c)). The hand chosen method using the BAD also categorizes from time to time atoms in the category "other", like in Fig 9 (b), but it is not able to capture the line defect.

The last method, the autoencoder applied to the BOO parameters, is also designed to detect automatically the different structures present within the system. However, when applying the method to our system, only two structures were detected: one associated with FCC and the other associated with HCP (the stacking fault). Also, only two structures were visible on the lower dimensional subspace over which the clustering was applied when following the method (see appendix Appendix B). We thus conclude that [START_REF] Ackland | Applications of local crystal structure measures in experiment and simulation[END_REF], (c) the CNA method [START_REF] Honeycutt | Molecular dynamics study of melting and freezing of small Lennard-Jones clusters[END_REF] from OVITO [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with OVITO -the Open Visualization Tool[END_REF], and (d) the combination of BOO and autoencoder as detailed in [START_REF] Boattini | Unsupervised learning for local structure detection in colloidal systems[END_REF]. We can remark on (a) two atoms belonging to the orange cluster from Fig 6 (c) within the dislocation lines. Atoms from this cluster were observed when two dislocations are close to each other. We can also remark on (b) one atom in red corresponding to the unassigned category on the hand-chosen criteria applied to BAD from [START_REF] Ackland | Applications of local crystal structure measures in experiment and simulation[END_REF].

the BOO parameters are not the most pertinent parameters to detect the structures within a plastically deformed crystal.

Figure 10: Snapshot of the system focusing on the interaction between two dislocations, one of them crossing the stacking fault of the other. These snapshots were obtained after filtering out the FCC structure. The atoms are colored depending on their surrounding structure obtained from (a) the combination of BAD and autoencoder, (b) hand-chosen criteria applied on BAD from [START_REF] Ackland | Applications of local crystal structure measures in experiment and simulation[END_REF], (c) the CNA method [START_REF] Honeycutt | Molecular dynamics study of melting and freezing of small Lennard-Jones clusters[END_REF] from OVITO [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with OVITO -the Open Visualization Tool[END_REF], and (d) the combination of BOO and autoencoder as detailed in [START_REF] Boattini | Unsupervised learning for local structure detection in colloidal systems[END_REF]. We can observe on (a) atoms belonging to the sky blue cluster from Fig 6 (c). Atoms from this cluster are observed after the interaction between a dislocation line and a staking fault.

Overall, these comparisons shows that our method is able to detect more structures than the other presented approaches. It is also able to capture the contours of the dislocation lines with a better precision than CNA as we can see on Fig 10.

Conclusion

In conclusion, we improved a method initially developed to detect automatically structures in colloidal materials. Through this method, we were able to finely detect different structure present in a Ni FCC crystal under plastic deformation. This algorithm uses the BAD parameters to describe the local environment of the atoms. Thanks to an autoencoder associated with a perturbation method, the most relevant BAD parameters were extracted. Then, by applying on these selected parameters two clustering methods, Kmeans and DBscan as well as the logistic regression classifier, different local structures have been successfully extracted. This procedure was able to detect the local structures present within the dislocation lines with a higher degree of precision compared to the commonly used CNA method. Our procedure is also more efficient than a method relying on hand-chosen criteria applied to BAD parameters or than combining autoencoder with BOO parameters.

Overall, this study shows that unsupervised learning applied to the BAD parameters is a promising approach to obtain more precise structure detection within crystalline materials. As the method detailed in this paper is based on parameters using hand-chosen ranges in its definition, it is not directly applicable to other crystalline system, especially alloys. While modified BAD parameters adapted for specific systems exist [START_REF] Amodeo | Atomistic Simulations of Compression Tests on Ni 3 Al Nanocubes[END_REF], developing a method to automatically find the optimal ranges to define the BAD parameters for any materials would allow to achieve an universal method for structural detection for crystalline materials. [START_REF] Boattini | Unsupervised learning for local structure detection in colloidal systems[END_REF] on our training data set. We calculate for each atom their BOO parameters using pyscal [START_REF] Menon | pyscal: A python module for structural analysis of atomic environments[END_REF]. We then obtained an optimal dimension of two at the bottleneck. Finally, we applied gaussian clustering at the bottleneck output detected two clusters: the yellow one corresponding to FCC structure and the blue one corresponding to HCP structure. As no other clusters can be seen, we observe that BOO is not able to separate the dislocation lines from the staking faults.

By applying this approach, we obtain an optimal bottleneck dimension equal to two. We then show on Fig B .12 a scatter plot showing for each atom their position in the space of the two bottleneck dimensions that we name {btneck1, btneck2}. By following the clustering method detailed in [START_REF] Boattini | Unsupervised learning for local structure detection in colloidal systems[END_REF], we extract two clusters: one in yellow corresponding to the FCC structure and one in blue corresponding to the HCP structure (in the stacking faults). We can see that no other clusters that could correspond to the structures within the dislocation lines are visible on Fig B .12. From this, we deduce that the BOO parameter is not pertinent to separate the dislocation lines from the other structures in a crystal under plastic deformation.

Figure 1 :

 1 Figure 1: (a) Snapshot made with Ovito (Stukowski, 2009) of the system used in this study. It consists of a cubic nanoparticle of nickel with edge length equals to 20nm. The two arrows represent the uniaxial compression performed with a flat indenter along the [100] direction. The stress-strain curve obtained from this compression is shown on (b). The vertical red line corresponds to occurence of the first plastic event.

Figure 2 :

 2 Figure 2: Schematic representation of autoencoder neural network architecture for dimensional reduction. The neural network is trained in order to reproduce the input at the output while passing through a bottleneck of lower dimensions. The encoder network aims at finding a low-dimension representation of the input. The decoder network is trained to reconstruct the original input from the lower-dimension bottleneck.

Figure 3 :

 3 Figure 3: Means square error E between the input and the output for different bottleneck dimensions. The error reaches a plateau for a bottleneck dimension of N = 3 which corresponds to the relevant number of BAD parameters sufficient to reconstruct the input data.

Figure 4 :

 4 Figure 4: Relative importance index of the BAD parameters obtained from the input perturbation method. The index is calculated by perturbating consecutively each input parameter and measure how much this perturbation affected the output. The parameters {χ 4 , χ 5 , χ 7 } are shown to be the most relevant ones to study the different structures in the system.

Figure 5 :

 5 Figure 5: The upper part of the figure (a) represents a Scatter plot of the BAD parameters projected on the most relevant parameters: {χ 4 , χ 5 , χ 7 }. The colors correspond to the substructures detected with K-means clustering. The lower part shows histograms showing the distribution of χ 4 (b), χ 5 (c), χ 7 (d). On the χ 4 histogram, we observe that most of the atoms are concentrated on two values, while on the χ 5 and χ 7 histograms, most atoms are concentrated on a single value. Thus, two combinations of {χ 4 , χ 5 , χ 7 } contains most of the atoms and corresponds in practice to the FCC and the stacking fault (HCP) structures.

  . While the Fig 6 (a) clearly shows three blocks, the two detected clusters seems to correspond to half of the central block. To understand what was detected by the K-means algorithm, we look on Fig 5 the distribution of χ 4 (b), χ 5 (c) and χ 7 (d).

  Fig 5 (a) that it totally fails in identifying two clusters, corresponding to different substructures, located further away from the central block. To detect them, we use another clustering method: the DBSCAN combined with the logistic regression classifier detailed in 2.3.3.2 and 2.3.3.3, respectively. We show its results on Fig 6 (b).

  Fig 6 (a), DBSCAN shown on Fig 6 (b)) as highlighted on Fig 6 (c). From this figure, we can see that we managed to detect a total of six different clusters: the dark blue and yellow ones from the central structures from the central block, containing the 99.8% of the atoms, corresponding to the main crystalline structures. The red and green cluster, corresponding to the two other aggregation of points clearly visible in the scatter plot of Fig 6 (c). Finally, the sky blue and the orange clusters, which only corresponds to few atoms, and detected only in few time steps.

  Fig 6 (c) will be colored in yellow on Fig 7. On Fig 7 (a), we can see that most of the atoms are contained in the yellow cluster which we can then associate with the FCC structure as our system is a Ni FCC nano-crystal.

Figure 6 :

 6 Figure 6: Scatter plots of the BAD parameters projected on the most relevant parameters: {χ 4 , χ 5 , χ 7 }. The colors correspond to the clustering methods used to detect different structures with: (a) K-means clustering, (b) DB scan clustering combined with logistic regression classification, and (c) a combination of the clusters obtained through (a) and (b). Each point corresponds to at least one atom, the distribution for each parameter can be seen on Fig 5.In total, six differents structures were detected through this combination of two clustering methods: FCC (yellow), HCP in the stacking fault (blue), upper and lower part of the dislocation lines (green and red), the interaction between dislocation line and stacking fault (skyblue) and the interaction between two dislocation lines (orange).

Figure 7 :

 7 Figure7: Snapshot of the system. The atoms are colored depending on their cluster as shown in Fig6 c). On panel (a), we observe that most of the atoms belongs to the yellow cluster which we associate with the FCC structure. On (b), we show the same snapshot after filtering out the FCC structure. We can observe that the combination of the two clustering methods is capable to separate the atoms linked with dislocation lines (green and red atoms) and those correlated with stacking faults (dark blue atoms).

  Fig 8, we only focus on one dislocation line and the associated stacking fault. From this figure, we can see that our method is able to capture with much more details the different structures present within a FCC crystal plastically deformed. To further emphasize the efficiency of our method, on Fig 9 we show a snapshot of the system with two dislocation lines about to interact and on Fig 10 we show the interaction between a dislocation line and a stacking fault.

Figure 8 :

 8 Figure8: Snapshot of the system. The atoms are colored depending on their surrounding structure obtained from (a) our combination of BAD and autoencoder, (b) hand-chosen criteria applied on BAD from[START_REF] Ackland | Applications of local crystal structure measures in experiment and simulation[END_REF], (c) the CNA method[START_REF] Honeycutt | Molecular dynamics study of melting and freezing of small Lennard-Jones clusters[END_REF] from OVITO[START_REF] Stukowski | Visualization and analysis of atomistic simulation data with OVITO -the Open Visualization Tool[END_REF], and (d) the combination of BOO and autoencoder as detailed in[START_REF] Boattini | Unsupervised learning for local structure detection in colloidal systems[END_REF]. Note that for the four structure detection methods shown here, we removed the surface atoms beforehand.

Figure 9 :

 9 Figure9: Snapshot of two dislocations about to interact. These snapshots were obtained after filtering out the FCC structure. The atoms are colored depending on their surrounding structure obtained from (a) the combination of BAD and autoencoder, (b) hand-chosen criteria applied on BAD from[START_REF] Ackland | Applications of local crystal structure measures in experiment and simulation[END_REF], (c) the CNA method[START_REF] Honeycutt | Molecular dynamics study of melting and freezing of small Lennard-Jones clusters[END_REF] from OVITO[START_REF] Stukowski | Visualization and analysis of atomistic simulation data with OVITO -the Open Visualization Tool[END_REF], and (d) the combination of BOO and autoencoder as detailed in[START_REF] Boattini | Unsupervised learning for local structure detection in colloidal systems[END_REF]. We can remark on (a) two atoms belonging to the orange cluster from Fig6 (c) within the dislocation lines. Atoms from this cluster were observed when two dislocations are close to each other. We can also remark on (b) one atom in red corresponding to the unassigned category on the hand-chosen criteria applied to BAD from[START_REF] Ackland | Applications of local crystal structure measures in experiment and simulation[END_REF].

Figure B. 12 :

 12 Figure B.12: Scatter plot showing the output of the 2-dimensional bottleneck obtained by applying the method from[START_REF] Boattini | Unsupervised learning for local structure detection in colloidal systems[END_REF] on our training data set. We calculate for each atom their BOO parameters using pyscal[START_REF] Menon | pyscal: A python module for structural analysis of atomic environments[END_REF]. We then obtained an optimal dimension of two at the bottleneck. Finally, we applied gaussian clustering at the bottleneck output detected two clusters: the yellow one corresponding to FCC structure and the blue one corresponding to HCP structure. As no other clusters can be seen, we observe that BOO is not able to separate the dislocation lines from the staking faults.

  and were optimised to differential crystal structures.

	BAD Parameter Minimum cos(θ jik ) Maximum cos(θ jik )
	χ 0	-1.0	-0.945
	χ 1	-0.945	-0.915
	χ 2	-0.915	-0.755
	χ 3	-0.755	-0.705
	χ 4	-0.705	-0.195
	χ 5	-0.195	0.195
	χ 6	0.195	0.245
	χ 7	0.245	0.795
	χ 8	0.795	1.0

Preprint submitted toDecember 22, 2022

Appendix A. Linking detected structures and dislocation lines Figure A.11: Snapshot of the system. On (a) the FCC atoms are filtered out to only keep the atoms in a stacking fault (in blue) and in a dislocation line (in green and in red). To understand the relative position between the green and the red atoms, we show on (b) the dislocation lines (in green) obtained from DXA [START_REF] Stukowski | Structure identification methods for atomistic simulations of crystalline materials[END_REF]. The green arrows on the dislocation lines corresponds orientation of the dislocation line. The black arrows corresponds to the burger vectors direction. The figures (c), (d) and (e) corresponds to zooms on the dislocation lines for easier visualization. We remark that when the dislocation line direction goes towards the trigonometric direction around the z-axis while the dislocation goes toward the center of the dislocation loop (as in (c) and (e)), the red atoms are above the green ones.

On Fig A.11 we show a comparison between the structures detected with our method (a) and the dislocation lines obtained from DXA (Stukowski, 2012) (b). On Fig A.11 (b) the black arrows correspond to the burger vector of each dislocation and the green arrows show the dislocation lines directions.

Appendix B. Autoencoder and clustering applied to BOO parameters

To compare more precisly the BOO and BAD parameters to detect local structure, we reproduced the method from [START_REF] Boattini | Unsupervised learning for local structure detection in colloidal systems[END_REF], which served as an inspiration for our method, and applied it to our data set. The method consist basically on calculation the BOO parameters for each atoms and training the autoencoder to find the optimal bottleneck dimension. Then, our method diverge from ours as they apply gaussian clustering directly to the output of the bottleneck.