Benoît Appolaire

Ludovic Héricher

Elisabeth Aeby-Gautier

Modelling of phase transformation kinetics in Ti alloys -Isothermal treatments

Keywords: titanium, modelling, diffusive transformations, nucleation & growth laws

We have developed a model which is able to describe the formation of high temperature microstructures in β metastable titanium alloys. It is based on the classical theory of nucleation, and on growth laws obtained by assuming that the processes are controlled by the diffusion of the alloying chemical species only. Two morphologies are described: grain boundary allotriomorphs, and Widmanstätten plates gathered in colonies growing from the allotriomorphs. Moreover, a statistical description of the nucleation conditions on β grain boundaries has been introduced in order to obtain the global kinetics of the transformation. The model predictions proved to be in good agreement with experimental results.

Introduction

The final microstructures of titanium alloys are designed through complex thermo-mechanical treatments. In order to predict the phase transformation kinetics and the amount of α phase with different morphologies, it is necessary to model the different successive transformations taking the temperature history into account. Moreover, the deformation occurring prior to the transformation has to be considered too [1].

In this study, we present a model where the transformation kinetics is described by nucleation and growth laws specific to the different morphologies, which is potentially able to deal with the effects of a prior deformation. We have applied this model to a multicomponent metastable β titanium alloy where one can observe at high temperatures (i.e. > 700 ˚C) the formation of grain boundaries allotriomorphs denoted α GB , and Widmanstätten plates colonies denoted α WGB , shown in Fig. 1. In section 2, we present the thermodynamic and kinetic data specific to the titanium alloys studied and necessary as inputs to the model. Then, the physical laws describing the nucleation and growth of the different morphologies observed are explained in section 3. The variation in the nucleation conditions due to the polycrystallinity of the material is taken into account when describing the global phase transformation kinetics in section 4. Finally, the kinetics of the global transformation as well as of some microstructural parameters are discussed in section 5.

Thermodynamic and kinetic data

The first requisite for modelling the phase transformation kinetics is a good description of the thermo-chemical equilibrium, i.e. the composition of the α and β phases and the corresponding tie-lines as a function of temperature. The equilibrium compositions have been obtained from samples heat treated at different temperatures for durations sufficient to achieve equilibrium in the β-Cez alloy [2], the composition of which is reported in Table 1. These experiments have shown that the composition of the α phase is quite independent of the temperature: thus, the α phase can be described as a stoichiometric compound. A great benefit can be drawn from this statement because the tie-lines are trivial. If this assumption appears very crude, it relies not only on our experimental measurements, but also on the following points. First, we have compared our measurements with the prediction of Saunders' database [START_REF] Saunders | Proceedings of the 8 th World Conference on Titanium[END_REF] input into Thermocalc software. This Calphad database encompasses the newest assessments of the thermodynamic parameters for titanium alloys [START_REF] Kaufman | [END_REF]5]. Large differences between both volume fraction of α and equilibrium compositions have been found, especially for temperatures lower than 800˚C, i.e. in the range of the isothermal treatments that we have studied. These differences led to strong deviations of the overall kinetics from the experimental IT diagram. Then, the use of the Calphad description to provide realistic kinetics would have required an optimization of the fitting parameters (e.g. the activity coefficients as well as the interaction parameters), a task which was outside the aims of our study. Second, our experimental evidence that the α phase composition does not vary much with temperature is not unique. This can be observed in the results of [6], where the variations in the compositions of V, Cr, Fe and Al (all elements present in the specific alloy studied) with temperature are very weak, and always lower than the reported experimental error bars. Finaly, a quick look at a few isothermal sections in some ternary titanium alloy diagrams shows that the curvature of the β phase field boundary is generally greater than that of the α phase field. Consequently, it is obvious that large variations in the equilibrium β compositions correspond to small changes in the α equilibrium compositions, due to the fan-shaped distribution of the tie-lines.

Assuming that the β phase is an ideal solution, the chemical equilibrium can be described with a solubility product [START_REF] Lupis | Chemical Thermodynamics of Materials[END_REF]:

i x β eq i x α i = K(T ) (1)
where x α i and x β eq i are the equilibrium molar fractions of the species i respectively in α and in β. Such a solubility product has been obtained from our measurements:

K(T ) = K 0 exp - Q R T ( 2
)
where R is the gas constant, K 0 = 0.819 and Q = 3500 J/mol. The assessed values give a very good description of the equilibrium α phase volume fraction as a function of temperature [START_REF] Héricher | [END_REF]. This solubility product is further used to calculate the driving force for nucleation and to determine the composition of the parent phase at the interface during the growth processes.

Because some diffusion data were lacking in the literature, we have performed experiments with diffusion couples made of pure Ti/β-Cez alloy, and determined the chemical diffusivities of the different species in the β temperature range [START_REF] Héricher | [END_REF].

3 Nucleation and growth laws

Transformation sequence

Some experimental studies have identified the main successive mechanisms governing the formation of α GB and α WGB at high temperatures [1,[START_REF] Bein | [END_REF][START_REF] Angelier | [END_REF]. First, α GB nucleates heterogeneously on β grain boundaries, and grows in directions normal to the grain boundaries. Second, as shown by systematic observations (Fig. 1b), some α WGB Widmanstätten plates appear along α GB allotriomorphs. Despite some observations in Ti-6.6 at.% Cr and Ti-8.6 at.% Mo alloys which have proposed sympathetic nucleation [START_REF] Kumar Menon | [END_REF], the specific appearance mechanisms of α WGB are still unclear. Third, we have clearly established that the α WGB lamellae grow inside the β grains, gathered in colonies where they have identical crystallographic orientations. These observations are in agreement with previous studies [12,13]. The β phase between the α lamellae is enriched in β stabilizer species. Furthermore, Laude [1] showed that the β phase is enriched in β stabilizer species in front of the colonies, over short distances which decrease with the transformation temperature (< 10 µm at 790˚C, < 7 µm at 750˚C). Consequently, the average composition of the β matrix outside the colonies can be assumed not to change during the transformation.

Based on these experimental observations, the model describes the transformation with the following sequence (Fig. 2

Nucleation of α GB

Despite the use of the pill box shape for describing the allotriomorph nucleation in some titanium based alloys [14], we have assumed that the α GB nuclei are spherical caps, in order to reduce the physical parameters involved in the model such as the interfacial energies, which are very difficult to assess experimentally.

The nucleation rate Ṅn of α GB (where N n is the number of nuclei and the dot denotes the time derivative) on the β grain boundaries is obtained using the classical theory of nucleation [15]:

Ṅn = N β Z exp -∆G * het kT exp (-τ /t) = Ṅns exp (-τ /t) (3)
where the different terms have their classical meaning and can be found in any standard textbook (see e.g. Ref. [START_REF] Christian | Theory of Transformations in Metals and Alloys -Part I[END_REF]); the only term specific to this study is the nucleation driving force expressed with the solubility product K (Eq.

2) as ∆G

c v = R T /V m ln i x β 0 i x α
i /K where V m is the molar volume and

x β 0 i the nominal molar fraction of element i in the alloy.

Growth of α GB

Based on the microstructural observations (Fig. 1a), the growth of α GB is assumed to be planar. Furthermore it is supposed to be controlled by the diffusion of the partitioning species. The width W of the α GB layers is calculated by solving the conservation equations for each partitioning species in the β matrix with the following boundary conditions (e.g. [17]). (i) For an homogeneous molar volume, the balances for each partitioning species i are:

x β * i -x α i Ẇ = J β * i (4)
where x β * i is the molar fraction of species i at the interface in β; and J β * i is the density of flux of species i at the interface in β. For an ideal solution J β i is simply given by the first Fick's law:

J β i = -D i (∂x β i /∂z). (ii)
The assumption of local equilibrium at the interface is described by the solubility product:

i x β * i x α i = K(T ) (5)
(iii) Moreover, the β matrix is considered as an infinite medium. This is valid as long as there is no change in the composition of the matrix far from the interface (equal to the nominal composition x β 0 i ), i.e. as long as the diffusion fields induced by the growth of competing α GB do not overlap. This is reasonable for α GB because the order of magnitude of the diffusion lengths in front of the α GB allotriomorphs is much smaller than the size of the β grains (a few hundred µm). The resulting equations to be solved are:

i (x β 0 i -x α i Ω i )/(1 -Ω i ) x α i = K(T ) Ω i = F η * 2 i = √ π η * i exp η * 2 i erfc (η * i ) (6)
where Ω i is the reduced supersaturation of the matrix in species i

Ω i = (x β * i - x β 0 i )/(x β * i -x α i ) and η * i is the width W scaled by the diffusion length of the species i η * i = W/(4D i t) 1/2 .

Appearance of α WGB

The appearance of α WGB is defined empirically because of the lack of knowledge concerning the mechanism involved. We have used the measurements of the final width of α GB W f as a function of temperature, realized by Laude [1]. A polynomial fit of the measurements in the temperature [600˚C-890˚C] gives:

W f = 7.98 • 10 -11 T 2 -9.59 • 10 -8 T + 2.93 • 10 -5 (7)
For temperatures lower than 600˚C, W f is assumed to be the same as for 600˚C (i.e. 0.488 µm) because of the lack of measurements.

Growth of α WGB

Because no data were available in the literature concerning the growth kinetics of α WGB in multicomponent β metastable titanium alloys, we have measured the average α WGB plates length on SEM images at different times of an isothermal treatment. A constant lengthening rate of L = 4.33•10 -8 ± 1.5 • 10 -8 m/s is found for a transformation at 790˚C. This value compares well with that found experimentally in Ti-6.6 at.% Cr by [14].

This observation supports the description of the α WGB diffusion controlled growth with a model accounting for the so-called point effect (e.g. see Ref.

[19] and references therein). We have assumed that the plates are parabolic cylinders which grow by a diffusion controlled process. As for α GB , the growth model relies on solving the conservation equations for each partitioning species in the β matrix. These are expressed in the stationary moving parabolic coordinate frame [START_REF] Ivantsov | translation of[END_REF]22] with the following boundary conditions: (i) Equation ( 4)balances at the interface for each partitioning species. (ii) The assumption of local equilibrium at the interface. It is also assumed that it is the diffusion at the very end of the plates, where the radius of curvature ρ is the smallest, which controls the lengthening kinetics. Under these conditions, the equilibrium is expressed as follows:

i x β * i x α i = K(T ) exp γ αβ V m ρ R T = K ρ (T ) (8)
(iii) The assumption that the β matrix is an infinite medium for reasons similar to those put forward for α GB . Indeed, the characteristic lengths of diffusion in front of the colonies are very small compared with the β grain sizes. The resulting equations are (see e.g. Refs. [START_REF] Ivantsov | translation of[END_REF]22]):

i (x β 0 i -x α i Ω i )/(1 -Ω i ) x α i = K ρ (T ) Ω i = F (P i ) (9)
where P i is the Péclet number of the plate edge

P i = ρ L/(2D i ).
Contrary to the case of the planar growth, there are two characteristic lengths for describing the plates, i.e. the length related to the growth rate L on one hand, and the width related to the radius of curvature ρ at the plate extremity on the other hand; thus, one more relation is necessary to solve the problem. When plotted against the radius of curvature, the velocity given by ( 9) reaches a maximum value. Zener [23] suggested the radius corresponding to this maximum velocity as the operating radius of the plates. Although empirical, this criterion proved to be quite successful for the Widmanstätten plates in other alloys (e.g. Refs. [24,25]). Equating to zero the partial derivative of ( 9) with respect to ρ gives the necessary closure equation.

Finally, because it is difficult to assess precisely the interfacial energies, the value of γ αβ has been fixed such that the calculated velocities fit to the measured ones over a large temperature range, with 790˚C as a guide: a value of 270 mJ/m 2 is obtained as the best fit. This value is in agreement with the nature of the semi-coherent interface which has been observed for that kind of morphology in binary titanium alloys [18,26].

Microstructural observations have shown that the α GB amount depends greatly on the β grain boundaries involved. This can be attributed to the difference in the α GB nucleation rates which are a function of γ ββ , i.e. the nature of the grain boundaries. To take into account such a scattering, it is necessary to introduce a statistical description of the transformation. For that purpose, we have considered one thousand tetrakaidecahedra (14-faced polyhedra) [27], representative of a polycrystalline β phase. Comparisons with different numbers of β grains have shown that one thousand is statistically sufficient, as well as computationally efficient for including the model in macroscopic mechanical computations of internal stresses [START_REF] Costa Teixeira | Proceedings of the International Conference on Thermal Process Modelling and Computer Simulation[END_REF]. All the tetrakaidecahedra have the same size, defined by the distance D between opposite faces. In the present study, the value of D is fixed to D = 200 µm, corresponding to the average diameter of the grains in the experimental samples obtained by image analysis. All the necessary geometric quantities, such as the volume and the surfaces are derived easily from D. This polycrystalline representation appears satisfactory when comparing the predicted specific surface S v = 17800 m -1 for D = 200 µm with the experimental one S v = 16900 ± 300 m -1 . A random wetting angle θ is allocated to each side of the tetrakaidecahedra faces within a Gaussian distribution (for simplicity) accounting for the scattering of the nature of the grain boundaries. The distribution parameters, i.e. the average value θ and the standard deviation δθ, are fitting parameters of the model. From preliminary calculations at several temperatures, these parameters have been set to the following values: θ = 19˚and δθ = 15˚. Thus we have γ ββ ( θ ∓ δθ) ∈ [500-540 mJ/m 2 ] all values included in the range [500-900 mJ/m 2 ] reported in Ref. [26] for titanium.

The implementation of the model is straightforward. All the faces (indexed j ∈ [1,14]) of all the grains (indexed k ∈ [1, 1000]) are treated sequentially at each time step in the following way. First, it is checked which transformation step is involved (Fig. 2). Then the appropriate microstructural variables are updated after having solved (3), ( 6) or ( 9): (i) The area covered by the α GB nuclei S jk n = π(r * jk sin θ jk )2 where the sum is performed over the total number N jk n of nuclei on the grain boundary j of grain k, i.e. N jk n = Ṅ jk n dt. The duration from the begining of the treatment necessary for the α GB nuclei to cover completely the face area is called the wetting time ∆t w . (ii) The width W jk of the α GB allotriomorphs layer. The corresponding transformed volume V jk GB = S jk ββ W jk can be deduced from it. (iii) Or the length L jk of the α WGB plates colony. The corresponding transformed volume V jk WGB = S jk ββ L jk is also computed. It must be emphazised that this transformed volume includes both α and β phases; the plates themselves, the width of which is supposed to be equal to the size of the face; and in between, the β phase enriched due to the thickening of these plates. For simplicity, it is assumed that both phases are in equilibrium inside the colonies. Then, the proportion of both phases is given by the lever rule. For each transformation step, it is checked whether it must be changed into the next one, or whether it must be stopped when the α WGB growth is concerned. The criteria for changing or stopping are: (i) Once 90% of the surface S jk ββ of the grain boundary j involved is covered by the nuclei (i.e. S jk n = S jk ββ ), the growth of the α GB allotriomorphs can start 2 . (ii) Once the width W jk of the α GB allotriomorphs layer has reached the experimental final width W f given by [START_REF] Lupis | Chemical Thermodynamics of Materials[END_REF], a new colony of α WGB can start to grow from the whole surface of the α GB layer (equal to the β grain boundary surface S jk

). (iii) The α WGB growth is stopped in two cases: first, if the total transformed volume j (V jk GB + V jk WGB ) is equal to the volume of the grain k; secondly, if the length L jk of the colony is greater than the grain size D minus the length of the colony growing from the opposite face of the same grain k, so as to prevent unphysical lengths of colony. (iv) The transformation is supposed to be finished when the transformation rate monitored by ξ is lower than 10 -6 s -1 , where ξ is the transformed volume fraction

V t /V tot = j,k (V jk GB + V jk WGB )/V tot .
Finally, we recall the microstructural input parameters used in the calculations in table 2. 5 Results

Isothermal transformation kinetics

The results of the calculations for an isothermal treatment at 790˚C are shown in Fig. 3. The transformed volume fraction ξ is plotted versus time in Fig. 3a. The corresponding microstructural parameters are plotted versus time in Fig. 3b: the average number of allotriomorphs layers per grain, the average number of colonies per grain, and the number of active grains where the transformation is progressing (i.e. where j V jk WGB > 0). The evolution of ξ can be split into two stages. During the first one, up to about 1000 s, a small increase is observed. It corresponds to the growth of α GB alone, as shown in Fig. 3b: the number of allotriomorph layers increases abruptly from the very begining and then slowly, whereas the number of α WGB colonies is null up to about 1000 s. Thus, ξ is nearly proportional to √ t, because the number of layers is quite constant, and because all the layers follow the same growth law W = C √ t with the same proportionality factor C. This stage ends abruptly when the width of a majority of the α GB layers reaches the empirical value W f . At the beginning of the second stage, there is a sharp increase in the number of α WGB colonies, which is simply the translation of the sharp increase of the number of allotriomorph layers at the begining of the treatment.

During the second stage, the growth of the α WGB colonies induces a rapid increase in the transformed volume fraction. Indeed, because the lengthening rate is time independent, the rate of transformation ξ is constant for a constant number of growing colonies. The slight convexity of ξ(t) ( ξ > 0) observed during the begining of this second stage (1000 s ≤ t ≤ 2500 s) is due to the slight increase in the number of colonies. For t > 2500 s, ξ slows down progressively ( ξ < 0) when the number of active grains starts to diminish. At the end of the transformation at 790˚C, the average number of colonies per grain is about five, lower than the maximum value imposed by the geometry of the tetrakaidecahedra.

In Fig. 4a, we have compared the global kinetics with the experimental in-situ measurements by electrical resistivity [1]. Because of the geometrical constraints, the value reached by the transformed volume fraction ξ may be less than unity at the end of the transformation (i.e. when ξ < 10 -6 s -1 ). Thus, we have plotted ξ ′ , i.e. ξ normalized by the transformed volume at the end of the transformation . The overall agreement is good. It must be emphasized that the duration of the first stage is close to the so-called incubation time of the order of 1000 s deduced from the experimental measurements. Usually, the experimental incubation time is attributed to the transient nucleation process (i.e. τ in (3)). But, at 790˚C, τ < 1 s for wetting angles θ < 45˚. Moreover, as stated above, during this first stage the degree of transformation is only a few per cent which is underneath the experimental uncertainty at the beginning of the measurements. Thus, it is reasonable to attribute the so-called experimental 'incubation time' to the growth of the α GB rather than to the nucleation process. Additional experiments were performed and have confirmed this result.

In Fig. 4b, we have compared the number of colonies vs. time predicted by the model and the measurements of Laude [1]. The predicted evolution is quicker and more sudden than the experimental one. Moreover, the calculated number of colonies at the end of the transformation is smaller than what has been measured, even when considering the experimental uncertainties of two colonies. To explain the time evolution of the α GB layers/α WGB colonies, we have plotted in Fig. 5 the wetting times ∆t w vs. the corresponding wetting angles θ at 790˚C. Two domains of θ are clearly separated by the angle θ c which is related to the duration of the transformation ∆t tot at 790˚C (dashed lines). For θ < θ c , ∆t w is small when compared to ∆t tot . In contrast, for θ ≥ θ c , ∆t w >> ∆t tot . Hence, as soon as a treatment starts, the faces with θ < θ c are immediately wetted. Then, the number of α GB layers increases suddenly. Only a few faces with θ ≈ θ c are wetted during the rest of the transformation; a slight increase in the number of α GB layers is then observed.

In Fig. 5, we have also plotted the distribution of the wetting angles for the Gaussian centered around θ = 19˚with δθ = 15˚. The proportion of the wetted faces, and so of the α WGB colonies at the end of the transformation corresponds to the shaded area divided by the total area under the Gaussian.

Influence of the temperature

The influence of the temperature is shown in figure 6. Concerning the first step, the higher the temperature, the longer is the first stage of transformation. It can also be noticed that the volume fraction of the α phase at the end of the first stage is the same whatever the temperature is, because the empirical criterion used for predicting the appearance of α WGB does not lead to large variations in the total volume of α GB [1]. Nonetheless, the slight increase of the final width W f with temperature is sufficient to lengthen the first stage as the temperature increases because W ∝ √ t. Similarly, during the second stage, the higher the temperature, the slower is the transformation rate. This is related to the number of colonies and to the lengthening rate L which both decrease with increasing temperatures, above 730˚C. Lastly, for the microstructural parameters chosen (Table 2), whatever the temperature is, the final volume fraction of α is slightly below the equilibrium one. This is a consequence of the incomplete fill-in of the grains by the colonies due to the geometrical constraints on their maximal lengths.

These calculated results are compared to experimental transformation kinetics obtained by electrical resistivity [1] for 750, 790 and 830˚C in Fig. 7a. The overall agreement is good. One can therefore perform additional calculations under isothermal conditions to construct the IT diagram of the β-Cez alloy (Fig. 7b). The calculated 10% and 90% of transformation are in good agreement with the experimental measurements in the high-temperature range. For lower temperatures (≤ 700˚C), the transformation kinetics obtained experimentally is clearly faster than what is predicted by the model. Indeed, the discrepancy is related to the appearance of a thin intragranular Widmanstätten structure (α WI ). Concerning the beginning time of transformation, the comparison is more complicated due to the uncertainty in the experimental measurements which ranges from typically ≈ 200 s at the highest temperatures to typically ≈ 30 s at the lower ones. In Fig. 7b, we have reported the calculated 1% and 2% of transformation. These curves show that, in the high-temperature range (> 700˚C), the first 3% of transformation span large time ranges in which the experimental beginning times of transformation are included. These large ranges can be explained by the change in the growth process between 1% and 2% of transformation: only α GB is predicted to grow at 1% of transformation, while α WGB is at 2%. Hence, the overall agreement between the model and the experiments can be considered as very satisfactory. Moreover, the model predicts the formation of α GB + α WGB morphologies at temperatures where a majority of α WI morphology appears. Experimental observations confirmed the occurence of these morphologies for the lower temperatures of 600 and 650˚C.

Besides the global kinetics, microstructural parameters can be predicted by the model. In Fig. 8 we have reported the average number of colonies per grain at the end of the transformation as a function of the temperature. The final number of colonies increases when the transformation temperature de- creases. Indeed, as the temperature decreases, the driving force for nucleation increases: the number of β grain boundaries which are wetted by α GB layers increases as well, as do the number of colonies which have grown from α GB . As shown in Fig. 8, these predictions agree qualitatively with the experimental observations. Nonetheless the predictions are lower than the measurements. In Fig. 9, we have compared the computed plate-edge radius with their experimental widths determined on SEM images with an overall uncertainty of about 150 nm, at different temperatures. The same trend is observed for both: they increase as the transformation temperature does. Furthermore, a factor of about 30 between the width and the radius is found whatever the temperature is.

Influence of parameters sensitive to a prior deformation

The thermal and mechanical treatments designed for these titanium alloys often start with a deformation step (e.g. forging). During this step, subgrains are formed that dynamic or metadynamic recrystallization can change into grains with small sizes, resulting in a grain refinement as well as a broadening of the grain boundary energy distribution [START_REF] Furuhara | [END_REF]. The microstructural parameters of the model which can take into account these effects are the distribution of the wetting angle θ, and the grain size D.

First, we have changed the width of the Gaussian distribution around the average value θ = 19˚. The time evolutions of ξ resulting from the calculations performed with different values of the standard deviation δθ are compared in Fig. 10. Two ranges of δθ are observed with opposite behaviours: for δθ < 10˚, when δθ increases the global kinetics accelerate mainly because the α WGB growth stage is accelerated; for δθ > 10˚, increasing δθ induces a small slow down of the α WGB growth stage. For explaining this two-fold behaviour, the growth kinetics of the different morphologies can be discarded because they do not change for a constant temperature. On the contrary, this behaviour is closely related to the wetting of the faces which controls the number of the growing colonies (Fig. 11). For very small δθ (e.g. δθ = 1˚in Fig. 10), many of the faces are close to the critical wetting angle θ c : the appearance of new colonies is progressive and so is the transformed volume fraction, leading to very low transformation kinetics. For moderate δθ ≤ 10˚, enlarging the Gaussian increases the number of faces which are rapidly wetted and thus increases the number of α WGB colonies (Fig. 11). This increase leads to a split of the global kinetics into the two identified growth stages as well as an acceleration of the α WGB growth stage. For large δθ > 10˚, a part of the Gaussian spans unphysical negative values of θ. These values are then discarded and new positive values following again a Gaussian distribution are allocated to the corresponding faces. The resulting distribution is the sum of two Gaussians, the area of the second one being the area of the part of the first one with the negative θ. This results in a skew distribution where the percentage of wetted faces is lower than for a symmetric Gaussian. The influence of the grain size is illustrated in Fig. 12. We have compared the global kinetics for different grain sizes D (100, 200 and 400 µm), keeping the other microstructural parameters equal to those of Table 2. A global acceleration is observed when the grain size diminishes. Moreover, the final transformed volume fraction ξ decreases when the grain size increases. The duration of the α GB growth is independent of the grain size because changing D does not change the final width W f of the α GB layers. However, the total volume fraction of α GB is larger for smaller grain sizes due to the larger specific surface of grain boundaries. Concerning the α WGB formation, changing the grain size has no effect on the growth rates Ẇ and L because the enrichment/depletion of the matrix is not taken into account. Moreover, in the present case, there is no effect on the number of colonies because the distribution of wetting times ∆t w is only a function of the temperature. Thus the rate ξ = Vt /V tot ∝ 1/D. However, for smaller grain sizes, it can be expected that the earlier growth of a few colonies may lead to an earlier complete filling of the grain volume, reducing the number of colonies per grain. Finally, the fact that the value of ξ reached at the end of the transformation depends on the grain size is quite unexpected. At first glance, this feature could be attributed to the adjustment of the matrix composition as transformation proceeds. But, even for large fractions of α phase, the average composition of the matrix does not vary much: first because the contribution of the allotriomorphs is very weak (at most 2% of the total volume is transformed at the grain boundaries); and secondly because the colonies of the Widmanstätten structure are composed of both phases with the respective equilibrium compositions. This last assumption is very reasonable because the diffusion length is of the order of the width of the α plates, i.e. at most a few µm (Fig. 9). Then, it can be safely assumed that the composition of the β phase in between the plates adjusts very quickly to reach the equilibrium, while the average matrix composition ahead of the colonies' advancing front is barely changed (in a similar way to what can happen in front of eutectoid colonies). Then the soft impingement of the colonies plays a very minor role in stopping the transformation, even if it could smooth the end of the transformation curves. For isothermal treatments, taking into account the change in the matrix composition is not relevant to explaining that the final ξ depends on the grain size. In fact, the reason why the model predicts such a dependency is the geometrical constraints that we imposed on the colonies. First, recall that we have considered constraints on their lengths as necessary for a good assessment of the mechanical properties. Indeed, the colonies could have been grown until the whole space was filled. Nonetheless, in that case, the colony lengths could reach a few grain diameters, which is not observed experimentally and which is problematic for assessing mechanical properties from the microstructure sizes. When the colonies are stopped for geometric reasons, untransformed 'dead zones' can remain if only a few faces of the grains have been wetted by the allotriomorph layers. The unwetted faces have very small nucleation potential due to their unfavorable energy. Some driving force can thus be still available for the transformation to proceed. If the stopping criterion on the rate of transformation is removed, the transformation proceeds at a rate controlled by nucleation on the unfavorable faces. This rate is so slow that the time for completion of the transformation is far greater than the duration of a typical experiment (several orders of magnitude larger). Because this dead-zone effect is accentuated as the grain size increases, the equilibrium fraction is all the more difficult to reach since the β grains are large. It must be noted that experimental results have shown the formation of fine α grains in such dead zones inherited from slow continuous cooling [30]. The fact that the colonies keep the same size fixed by the tetrakaidecahedron grain geometry is the second geometrical constraint which contributes to the existence of dead zones. The next step to improve the modelling of such transformation sequences would be to take into account nucleation at the edges and corners. This has not been done yet because more experimental evidence for the relationship between the shapes of the colonies and the places where they nucleate are needed. In conclusion, it must be stressed that the previous results relative to both parameters δθ and D are still in qualitative agreement with the experimental in-situ measurements of the global transformation kinetics performed by Laude [1] after a prior deformation in the β phase field. Indeed, these exper-iments have shown that the kinetics are accelerated by a prior deformation. This agreement seems to support the relevance of our model. More work is currently in progress to assess quantitatively the change in the parameters of the model resulting from a prior deformation.

Conclusions

We have developed a model able to predict the kinetics of α GB and α WGB morphologies during isothermal treatments in some titanium alloys. It is based on nucleation and growth laws, and the necessary data are the thermodynamic description of the phases, the diffusion coefficients, and some structural parameters such as the interface energy, the grain boundary energies distribution and the average grain size. This model has been applied to the β-Cez alloy. The calculated overall transformation kinetics are in good agreement with experimental measurements. Moreover, the evolution of some microstructural parameters, which are important for the mechanical properties such as the mean number of colonies per grain, follows the same trends as observed experimentally. Finally, provided that the thermodynamic description of the phases is available, this model can easily be used for other titanium alloys where similar sequences of morphologies are observed. Furthermore, it can take into account the effect of a prior deformation of the β phase which is known to modify the nucleation conditions.

Fig. 1 .

 1 Fig. 1. SEM images of the α phase morphologies observed at 790˚C in β-Cez alloy: (a) α GB and (b) α WGB colonies growing from α GB .

Fig. 2 .

 2 Fig. 2. Sequence of the phase transformation modelling for the high temperature morphologies: (a) nucleation of α GB , (b) growth of α GB , (c) appearance of α WGB , (d) growth of α WGB .

Fig. 3 .

 3 Fig. 3. Isothermal treatment at 790˚C: (a) transformed volume fraction ξ vs. time; (b) average numbers of α GB allotriomorph layers and of α WGB colonies per grain, and number of active grains vs. time.

Fig. 4 .

 4 Fig. 4. Comparison with experimental results at 790˚C [1]: (a) degree of transformation ξ ′ vs. time; (b) average numbers of α WGB colonies per vs. time. Continuous lines are for the calculations and dots for the measurements.

Fig. 5 .

 5 Fig. 5. Wetting time and density of faces per angle vs. wetting angle at 790˚C.

Fig. 6 .

 6 Fig. 6. α volume fraction vs. time for different temperatures. The dashed lines are the equilibrium α phase volume fractions.

Fig. 7 .

 7 Fig. 7. Comparison with experiments: (a) global kinetics for 750, 790 and 830˚C; (b) IT diagram of β-Cez in the high temperatures range. The lines are for the calculations and the dots for the experiments. For clarity, the error bars have been omitted for the experimental 10% of transformation.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Average final number of colony vs. temperature.

Fig. 10 .

 10 Fig. 10. Time evolution of ξ for different standard deviations δθ at 790˚C. The arrow follows the increasing δθ

Fig. 11 .

 11 Fig. 11. Number of α WGB colonies vs. the standard deviations δθ.

Fig. 12 .

 12 Fig. 12. Degree of transformation ξ vs. time for different grain sizes D.

Table 1

 1 Composition of the so-called β-Cez titanium alloy.

	Species	Ti	Al Mo Cr Fe Zr Sn
	at. %	balance 8.8 2.	2. 0.8 2. 0.8

  WGB at the growing α GB front; (d) growth of α WGB from α GB .

	): (a) nucleation of α GB on a β
	grain boundary; (b) growth of α GB on one side of a β grain boundary; (c)
	appearance of α

Table 2

 2 Microstructural input parameters.

	Values	200	270	19 15

Parameters D (µm) γ αβ (mJ/m 2 ) θ δθ

Corresponding author. Fax: (33) (0)3 83 58 40 56. e-mail: appolair@mines.inplnancy.fr

Neglecting the incubation time, the integration of (3) gives a infinite time for exhausting all the nucleation sites. This result is unrealistic because the nuclei are able to grow along the grain boundaries, decreasing the number of nucleation sites. The 90% criterion is a first way to account for this process. Some additional work is under process to further analyse the wetting of the grain boundaries by allotriomorphs.


L Hervo

J.-M Senoner

A Biancherin

Large-Eddy Simulation of Kerosene Spray Ignition in a Simplified Aeronautic Combustor

The current work presents the Large Eddy Simulation (LES) of a kerosene spray ignition phase in a simplified aeronautical combustor for which detailed experimental data are available. The carrier phase is simulated using an unstructured multi-species compressible Navier-Stokes solver while the dispersed liquid phase is modeled with a Lagrangian approach. An energy deposition model neglecting the presence of a plasma phase in the very first instants of the energy deposition process, a reduced kinetic scheme and a simplified spray injection model are combined to achieve both a reasonable computational expense and a satisfactory overall accuracy. Following a brief description of the validation of these models, non reactive gaseous and twophase flow LES's of the target combustor are performed. Excellent agreement with experiments is observed for the non reactive gaseous simulations. The dispersed phase velocity fields are also well reproduced while discrepancies appear for the spatial size distribution of the particles. Finally, numerical snapshots of a successful ignition phase are shown and discussed.

Introduction

In order to reduce the consumption and the pollutant emissions of their engines, aircraft manufacturers are testing new combustor designs that may be operated safely in a lean regime. However, lean premixed combustors may be difficult to ignite, especially at high altitude conditions. Today, numerous experimental tests are required to ensure ignition capability in the whole operating domain. In this context, the availability of accurate numerical tools could provide better physical insight into ignition phenomena, shorten the design cycles of new combustors and reduce the necessity of expensive test campaigns. However, simulating the ignition phase of an aeronautical combustor is a difficult task as this event is a complex process involving a variety of physical phenomena. Five main phases may be distinguished during an ignition sequence [22]:

-The first phase is characterized by the energy deposition process. From a technological point of view, ignition is still mostly initiated with spark plugs in modern aircraft engines, despite important research efforts on alternative techniques such as laser focusing, micro-waves, etc... Regardless of the device, the energy deposition is generally very localized. For aeronautical applications, local energy concentration levels are such that temperature levels exceed several tens of thousand Kelvins, leading to a local ionization of the atoms in the gas and the generation of a plasma phase. Energy deposition is also characterized by a very fast and significant pressure increase, typically of the order of several bars. -As the pressure suddenly increases, a shock wave contributing to the expansion of the plasma generally forms. As a consequence of this expansion, the maximum temperature and pressure levels start to decrease. -As temperature levels decrease to several thousands of Kelvins, a kernel of hot gases appears. Provided the kernel's chemical composition is favourable, it may ignite and form a flame kernel. -The evolution of the flame kernel is then dictated by its propagation and its convection within the aeronautical combustion chamber. Furthermore, the flame kernel becomes wrinkled, stretched and strained as it evolves in a turbulent flow field. For several reasons, it may be quenched during this phase ("failed ignition") or continue to expand within the chamber until a stabilized combustion regime is attained ("successful ignition"). Therefore, ignition in turbulent flow configurations may only be defined via local probabilities. -Additionally, real combustors are made of several burners and ignition is generally not triggered in all of them. Thus, flame propagation to all burners must also be ensured.

For single burner configurations, the main quantity of interest is the local ignition probability as a function of the operating conditions. In a first approximation, this quantity may be estimated from local properties such as the equivalence ratio, the mean and fluctuating velocity fields, etc... Models based on Reynols Averaged Navier Stokes (RANS) or time averaged Large Eddy Simulations (LES) as well as instantaneous LES flow fields or a combination of both may be found in the literature [START_REF] Eyssartier | Using LES to predict ignition sequences and ignition probability of turbulent two-phase flames[END_REF][START_REF] Neophytou | Spark ignition of turbulent recirculating non-premixed gas and spray flames: A model for predicting ignition probability[END_REF]. Due to the continuous increase of computational resources, LES's of the ignition sequence of entire annular combustors are becoming feasible [START_REF] Boileau | LES of an ignition sequence in a gas turbine engine[END_REF]. Similarly, ignition probabilities may be directly evaluated from LES's [12].

From the above phenomenological description, it follows that these simulations would need to explicitly account for the presence of a plasma in the first instants following the energy deposition process. However, this requires coupling a Maxwell and a compressible Navier-Stokes solver [START_REF] Thiele | Numerical simulation of spark ignition including ionization[END_REF] and is rarely done in practice. Instead, the influence of the plasma phase on the properties of the kernel of hot gases is often modelled. Duclos and Colin [START_REF] Duclos | Arc and kernel tracking ignition model for 3D spark ignition engine calculations[END_REF] proposed the Reynolds Averaged Navier-Sokes (RANS) AKTIM model, where the temporal evolutions of tension and current intensity within the discharge circuit are used to estimate the amount of energy supplied to the spark plug. The spark is then represented with Lagrangian tracers in order to account for the interaction with the surrounding gas. These tracers are initially carrying a vanishing amount of burnt gases, a given amount of fresh gases and an excess of energy. Whenever the excess of energy lies above a fixed threshold, the mass of burnt gases becomes equal to the mass of fresh gases and a flame kernel is initiated in the flow field. A very similar model relying on the use of Lagrangian tracers until the formation of the initial flame kernel named Spark Channel Ignition Monitoring Model (SparkCIMM) was proposed by Dahms et al. [START_REF] Dahms | Understanding ignition processes in spray-guided gasoline engines using high-speed imaging and the extended spark-ignition model SparkCIMM. Part A: Spark channel processes and the turbulent flame front propagation[END_REF]. Once this flame kernel exceeds a characteristic size allowing its explicit resolution on the numerical grid, the SparkCIMM model is deactivated and tracked with a combustion model based on the G-equation. The AKTIM model was later adapted from the RANS to the Large Eddy Simulation (LES) context by Colin and Truffin [START_REF] Colin | A spark ignition model for large eddy simulation based on an FSD transport equation (ISSIM-LES)[END_REF] and coined Imposed Stretch Spark Ignition (ISSIM) LES model. It relies on the addition of an energy source term to the ECFM combustion model [START_REF] Richard | Development of LES models based on the flame surface density approach for ignition and combustion in SI engines[END_REF]. Curvature effects occurring at the subgrid scale rely on an analytical model as long as the kernel is underresolved on the numerical grid. A simpler level of modeling may be achieved by imposing a temporally and spatially varying source term to the energy equation, as done in the Energy Deposition (ED) model proposed by Enaux [START_REF] Enaux | Simulation aux Grandes Echelles d'un moteur à allumage commandéévaluations des variabilités cycliques[END_REF] and Lacaze [START_REF] Lacaze | Large eddy simulation of spark ignition in a turbulent methane jet[END_REF]. In order to limit the influence of the subgrid scales in the first instants of the flame kernel formation, the authors combined their ED model with a high level of refinement in the area of energy deposition. Despite its simplicity, this approach was able to quantitatively predict the dynamics of a methane jet ignition [START_REF] Lacaze | Large eddy simulation of spark ignition in a turbulent methane jet[END_REF].

When the fuel is supplied to the kernel of hot gases as a spray, droplet flame interactions need to be additionally taken into account. The characteristic fuel droplet evaporation time scale will generally be much larger than the characteristic chemical time scale of combustion, leading to a decrease in flame propagation speeds compared to the perfectly premixed case [START_REF] Neophytou | Simulations of laminar flame propagation in droplet mists[END_REF]. Moreover, spray evaporation acts as a significant heat sink and may thus favour flame extinction phenomena. More details on spray ignition may be found in the reviews by Aggarwal [START_REF] Aggarwal | A review of spray ignition phenomena: Present status and future research[END_REF] and Mastorakos [START_REF] Mastorakos | Forced ignition of turbulent spray flames[END_REF].

The present work aims at assessing the capabilities of an energy deposition model, reduced chemical kinetics and a global spray injection model to simulate the ignition phase of a simplified aeronautical combustor in the presence of a dispersed phase. To the authors' knowledge, numerical simulations of the entire ignition phase have not yet been published for this target configuration despite previous investigations [START_REF] Bruyat | Influence de l'évaporation de gouttes multicomposant sur la combustion et des effets diphasiques sur l'allumage d'un foyer aéronautique[END_REF][START_REF] Eyssartier | Using LES to predict ignition sequences and ignition probability of turbulent two-phase flames[END_REF].

The current article is organized as follows. First, the simplified methodology used to model the energy deposition phase of the ignition process is described. Second, the target combustor is presented. The main characteristics of the Navier-Stokes solver used to simulate the non-reacting gaseous flow field within the target configuration are then detailed before presenting numerical results. Third, simulations of the dispersed two-phase flow within the target configuration are shown. Then, the simulation of a successful ignition phase is analysed. The main findings are summarized in the conclusion, followed by possible outlooks.

Ignition modeling

Energy deposition model

In order to approximately reproduce the phenomena leading to the formation of a kernel of hot gases in the vicinity of the electrodes of the spark plug, the Energy Deposition (ED) of Enaux [START_REF] Enaux | Simulation aux Grandes Echelles d'un moteur à allumage commandéévaluations des variabilités cycliques[END_REF] and Lacaze [START_REF] Lacaze | Large eddy simulation of spark ignition in a turbulent methane jet[END_REF] is used. This model simply consists in adding an unsteady source term to the energy equation. This source term is assumed to follow a Gaussian curve in time and space [START_REF] Enaux | Simulation aux Grandes Echelles d'un moteur à allumage commandéévaluations des variabilités cycliques[END_REF].

Q = i 4π 2 σ 3 S σ t e -1 2 r-r 0 σ S 2 e -1 2 ( t-t 0 σ t ) 2 (1)
with:

σ S = ∆ S
2 ln(10 4 ) and

σ t = ∆ t 2 ln(10 4 ) ( 2
)
where Q is the energy source term, i is the total energy supplied to the gas by the spark plug, ∆ S and ∆ t are the length and duration where 99% of the energy is deposed thanks to the normalization constant 2 ln(10 4 ), r is the distance from the center of the source term r 0 , t is the time and t 0 is the instant where the power given by the spark is maximum. It was verified that the implemented numerical procedure allowed to exactly control the amount of energy supplied to the flow field over time (not shown).

Thermodynamics and transport properties

The energy supplied to the flow by the aforementioned source term leads to significant local temperature increase. As a consequence, the range of validity of most thermophysical and transport properties may be exceeded. In order to extend this range of validity, the thermodynamic properties (heat capacity, enthalpy and entropy) of the major components N 2 and O 2 are described by polynomial expressions provided by McBride et al. [START_REF] Mcbride | NASA Glenn coefficients for calculating thermodynamic properties of individual species[END_REF] and valid up to 20 000 K in the present work. However, it is important to note that these polynomial expressions only account for temperature dependencies related to the molecular degrees of freedom although the aforementioned properties also depend on the chemical composition of the mixture. Accounting for composition changes due to endothermic ionization and dissociation processes would require dedicated kinetic mechanisms characterizing these phenomena and was considered out of the scope of the present work. Fig. 1 displays the comparison between the modified heat capacity curve at constant pressure for nitrogen N 2 [START_REF] Mcbride | NASA Glenn coefficients for calculating thermodynamic properties of individual species[END_REF] used in the present numerical simulations and the reference curves available in the CEDRE software (more details on the CEDRE software are provided in section 4.1). Moreover, a comparison between the theoretical polynomial curves and the actual values retrieved from numerical simulations with imposed temperature variations is also shown for validation purposes. Unfortunately, the authors could not find any data to describe the thermodynamic properties of combustion products such as carbon monoxide CO, carbon dioxide CO 2 or the various fuel species at temperatures above 5000 K. The latter were therefore kept constant above this threshold value.

The transport properties of the species are described by laws commonly used and no specific modifications were applied for higher temperatures. In particular, the viscosity of the different species follows a Sutherland law while the thermal conductivity and the diffusivity of the species are calculated with the assumption of a constant Prandtl and Schmidt number. It is important to note that the approximations made for transport properties must remain consistent with those used to derive the reduced kinetic schemes to ensure a correct laminar flame speed. However, it should be possible to combine these transport properties with different laws valid at higher temperatures and this aspect will be investigated in future work.

Validation case in laminar premixed conditions

In order to assess the capabilities of the aforementioned modeling ingredients on an elementary test case, the propagation of a laminar flame kernel in a propane-air mixture is simulated. The test case reproduces the conditions of an experiment performed by Renou and Boukhalfa [START_REF] Renou | An experimental study of freely propagating premixed flames at various Lewis numbers[END_REF], and was simulated by Enaux [START_REF] Enaux | Simulation aux Grandes Echelles d'un moteur à allumage commandéévaluations des variabilités cycliques[END_REF] with the AVBP solver. Chemical reactions are modeled with a two-step scheme reproducing correct laminar flame speeds at ambient conditions [START_REF] Enaux | Simulation aux Grandes Echelles d'un moteur à allumage commandéévaluations des variabilités cycliques[END_REF]. The parameters related to the energy deposition model, see eq. 1, are strictly identical to those of Enaux [START_REF] Enaux | Simulation aux Grandes Echelles d'un moteur à allumage commandéévaluations des variabilités cycliques[END_REF] and summarized in table 1.

∆ S 2 mm ∆t 120 µs i 2.5 mJ
Table 1: Parameters of the energy deposition model in the numerical simulations of the laminar flame kernel propagation experiment of Renou and Boukhalfa [START_REF] Renou | An experimental study of freely propagating premixed flames at various Lewis numbers[END_REF] The evolution of the flame radius over time is compared to experiment in fig. 2 where the results of Enaux [START_REF] Enaux | Simulation aux Grandes Echelles d'un moteur à allumage commandéévaluations des variabilités cycliques[END_REF] are also reported. The flame radius is defined as follows:

R = 12 π V T 1 - Y mean f Y 0 f 1/3 (3)
where V T , Y mean f and Y 0 f respectively denote the total volume of the computational domain, the mean and initial fuel mass fraction in this domain. As may be seen on fig. 2, the agreement between both simulations and the experiment is very good for the temporal evolution of the flame radius. An equally good agreement is observed for the flame speed, as may be seen on fig. 3. The flame propagation speed is defined as the temporal derivative of the radius according to eq. 3:

S c = ρ 1 ρ 0 dR dt (4)
with the subscripts 0 and 1 respectively denoting properties in the fresh and burnt gases. From fig. 3, it appears that the flame propagation is extremely fast until approximately t = 0.3 ms after which it starts decelerating at the expected laminar flame speed of approximately 0.2 m s -1 at these conditions. The very fast fuel consumption in the first instants may either be due to auto-ignition in the expanding kernel of hot gases or the high stretching of the flame kernel. Fig. 3: Evolution of the flame propagation speed over time for the ignition of a laminar propane-air mixture : experimental data from Renou and Boukhalfa [START_REF] Renou | An experimental study of freely propagating premixed flames at various Lewis numbers[END_REF], AVBP simulation from Enaux [START_REF] Enaux | Simulation aux Grandes Echelles d'un moteur à allumage commandéévaluations des variabilités cycliques[END_REF] and CEDRE simulation from the present work with with two values of deposited energy.

At this stage, it is important to stress that using such simplified energy deposition model requires an accurate estimate of the energy amount supplied to the gas by the spark. Unfortunately, uncertainties relative to the energy deposition typically amount to 30% or more [22,[START_REF] Teets | Calorimetry of ignition sparks[END_REF]. The amount of energy supplied to the flow is expected to determine flame propagation in the very first instants. To verify this assumption the simulation has been run for two values of total amount of energy supplied to the flow, 2.5 mJ and 1.25 mJ (i.e. 50%). The resulting shift is relatively minor in the laminar case, as may be seen from fig. 2. However, in a turbulent flow field and a realistic combustor geometry, the energy supplied to the gas represents a fundamental parameter to characterize the success or failure of an ignition attempt. Therefore, these uncertainties on the amount of energy supplied to the flow represent an important obstacle to the accurate prediction of ignition probabilities. While the presented methodology could be relatively accurate in the characterization of flame propagation, its use for the prediction of ignition probabilities seems questionable due to the important simplifications and uncertainties related to the characteristics of the kernel of hot gases resulting from the energy deposition process. For this reason, the present work will mainly focus on flame propagation dynamics during a successful ignition phase of the target configuration.

The Mercato configuration

The Mercato test rig was designed at ONERA to investigate two-phase combustion in a configuration representative of a realistic helicopter engine. This experimental facility can be operated at pressures between 0.5 bar and 1 bar and at temperatures ranging from 233 K to 473 K for the gas, allowing to reproduce high altitude conditions. The fuel is liquid kerosene and is injected at temperatures ranging from 233 K to ambient temperatures. The combustion chamber has a rectangular shape (130 mm x 130 mm) and it is equipped with Quartz windows for optical access. The air is injected into the plenum via a channel before passing through the swirler, where a rotational motion is imposed to the air. The fuel is injected at the chamber inlet via a Simplex pressure-swirl atomizer. Finally, a spark plug mounted on the wall of the combustion chamber is used for ignition. A view of the experimental setup is provided on fig. 4, while a schematic of the test rig is displayed on fig. 5. Although simplified, the Mercato configuration displays a swirler, an injection system and a spark plug which are representative of realistic aeronautical engines.

Multiple experimental campaigns were conducted for this test rig and data characterizing the purely gaseous flow and the two-phase flow in both nonreactive and reactive conditions are available [START_REF] Lecourt | TIMECOP-AE WP2 D2.2.1c -Injection system two-phase flow characterisation (LDA-PDA)[END_REF][START_REF] Rosa | Phénomènes d'allumage d'un foyer de turbomachine en conditions de haute altitude[END_REF]. Furthermore, ignition probabilities and trajectories of flame kernels in the combustion chamber were measured [START_REF] Lang | Statistical evaluation of ignition phenomena in turbojet engines[END_REF]21]. In our study, the so called reference point (see table 2) was chosen because it includes various data on both non-reactive and reactive gaseous flow and on the dispersed phase for the non-reactive case in addition to data characterizing ignition sequences. An implicit second-order Runge-Kutta time integration method (RKI2) and a second-order accurate MUSCL approach in combination with a Roe scheme are used for spatial discretization. The computational domain comprises the air inlet channel, the plenum, then the combustion chamber and the suction spout, see fig. 6. In order to partially compensate for the numerical dissipation of the Roe scheme, the numerical grid was refined with respect to a reference grid used to perform LES's of the same configuration with the AVBP solver [START_REF] Eyssartier | Study and modelisation of stationnary and transient two-phase flow combustion[END_REF].

The zones of highest refinement are located inside the swirler and the first third of the combustion chamber, see fig. for respectively mean and root-mean square (r.m.s) axial and tangential air velocity profiles. Equally good agreement is observed for the radial velocity profiles but not shown for concision. The comparison is performed at three different planes respectively located at axial distances z = 6, z = 26 and z = 56 mm from the chamber inlet plane, see fig. 6. The present results are in very good agreement with experimental data as well as numerical results obtained with the AVBP solver. In particular, the opening of the jet and the magnitude of the swirling motion are very well captured. Despite some lack of convergence, velocity fluctuation levels are also in very good agreement with experiments up to the third measurement plane z = 56 mm. The third measurement plane is where the energy deposition will be performed for the ignition simulations Fig. 6: Visualization of global mesh and refinements in a midplane cut along the x-axis of the computational domain and a zoom on the midplane cut through the swirler and the first third of the combustion chamber. The location of the experimental measurement planes, i.e. at z=6 mm, z=26 mm and z=56 mm is highlighted on the zoomed view.

and the present agreement indicates that fluctuations levels of the gaseous phase will be well reproduced around this region.

Two-phase flow

Injection procedure

The dispersed phase is described with the Lagrangian solver of the CEDRE platform named SPARTE [START_REF] Murrone | Numerical modeling of dispersed two-phase flows[END_REF]. In the Lagrangian framework, particles are treated as pointwise inclusions without explicitly resolving their surrounding flow fields. The underlying spray density function is then approximated through Dirac delta functions centred on each inclusion [START_REF] Williams | Spray combustion and atomization[END_REF]. In the present application, the particles result from the atomization of a swirled liquid sheet. First attempts to couple the dispersed solvers of CEDRE with a diffuse interface tracking technique to model liquid atomization are currently being made [START_REF] Zuzio | Multi-solver LES simulation of the atomization of a cross-flow liquid jet in a channel[END_REF], but they are not yet mature enough to be used for the ignition simulations performed in the present work. Therefore, particles are directly injected at the atomizer orifice using semi-empirical relations to specify their velocities while relying on experimental data to determine their size distribution [START_REF] Sanjosé | Fuel injection model for Euler-Euler and Euler-Lagrange large-eddy simulations of an evaporating spray inside an aeronautical combustor[END_REF]. The velocity profile used for particle injection is then made of two components, a constant axial velocity deduced from liquid mass conservation at the atomizer nozzle and an azimuthal velocity component resulting from the swirling motion imposed to the liquid within the atomizer. This swirling motion generates an air core along the axis of the atomizer such that liquid is only present between the radius of the air core R a and the orifice radius R 0 [START_REF] Lefebvre | Étude expérimentale et numérique de l'allumage des turboréacteurs en conditions de haute altitude[END_REF]. The resulting velocity profile then writes:

u l,z = ṁl ρ l πR 2 0 (1 -X) (5)
u l,φ = ṁl ρ l A p R s r ( 6
)
α 0 l (r) = 0 if r ∈ [0, R a ] 1 if r ∈ [R a , R 0 ] ( 7
)
where R 0 s = (R 0 + R a )/2 denotes the center radial position of the liquid at the atomizer orifice. The two remaining unknowns are the area ratio of the air core Table 4: Models used to describe the dispersed phase to the nozzle X and the area of the tangential inlet ports inducing the swirling motion A p . Both are deduced using empirical correlations of the literature [START_REF] Lefebvre | Étude expérimentale et numérique de l'allumage des turboréacteurs en conditions de haute altitude[END_REF] based on the orifice radius R 0 = 0.5 mm, the liquid mass flow rate ṁl = 2.25 g s -1 and the spray half angle θ 0 s = 40 • . The reader is referred to Sanjose et al. [START_REF] Sanjosé | Fuel injection model for Euler-Euler and Euler-Lagrange large-eddy simulations of an evaporating spray inside an aeronautical combustor[END_REF] for more details on the injection procedure. The experimental particle size distribution is approximated via a log-normal law:

f n (r) = 1 √ 2πσ ln r exp - ln(r/r ln ) 2 √ 2σ ln (8)
with r ln and σ ln respectively the droplet radius and the distribution variance.

For the present simulations, the values r ln = 5.76•10 -6 m and σ ln = 0.564 were used. This simplified injection procedure was validated separately by injecting particles in a quiescent environment and verifying the imposed velocity profiles as well as the resulting size distribution (not shown).

Numerical setup

To evaluate the particle drag at higher particle Reynolds numbers, the correlation of Schiller-Naumann [START_REF] Schiller | A drag coefficient correlation[END_REF] is applied. A standard Spalding evaporation model [START_REF] Spalding | A standard formulation of the steady convective mass transfer problem[END_REF] is used to evaluate heat and mass exchanges between the particles and the surrounding carrier phase. Convective effects are accounted for via the correlations of Ranz and Marshall [START_REF] Ranz | Evaporation from drops[END_REF]. In addition, the thickening of the boundary layers at the drop's surfaces due to the evaporative mass flux is accounted for using the relations of Abramzon and Sirignano [START_REF] Abramzon | Droplet vaporisation model for spray combustion calculations[END_REF]. Two-way coupling is enabled via the standard point force approximation [START_REF] Maxey | Localized force representations for particles sedimenting in Stokes flow[END_REF]. The time advancement of the particles is based on a semi-analytical first order scheme. The interpolation of gaseous properties at the particle locations relies on the gradients calculated by the carrier phase solver CHARME, yielding secondorder spatial accuracy. Finally, the particles are located on the computational grid using the algorithm of Haselbacher et al. [START_REF] Haselbacher | An efficient and robust particle-localization algorithm for unstructured grids[END_REF]. The main numerical parameters for the dispersed phase are summarized in table 4.

Results

The numerical velocity and diameter fields averaged over approximately 350 ms for the dispersed phase are compared to experimental data from Linassier [21] and simulation results obtained by Eyssartier with the AVBP solver [START_REF] Eyssartier | Study and modelisation of stationnary and transient two-phase flow combustion[END_REF] on fig. 9. It is important to note that the particle velocity data was obtained via a mass average. Also note that due to particle impingement on the windows of the chamber, a liquid film formed there and experimental data could only be acquired at a single plane located at z = 6mm from the chamber inlet plane. Despite a slight lack of convergence, numerical velocity profiles appear in good agreement with experimental data for the present work. As may be seen on fig. 10, numerical results are less satisfactory for the diameter distribution. Only the qualitative behaviour of an increase of the particle diameter towards the axis of the chamber is reproduced while quantitative levels are mispredicted. The overly simplified injection procedure could partly explain these differences. However, Senoner [START_REF] Senoner | Simulations aux grandes échelles de l'écoulement diphasique dans un brûleur aéronautique par une approche Euler-Lagrange[END_REF] obtained good agreement between numerical particle diameter profiles and experiment using the same simplified injection procedure for the hot operating point (see table 2). Also, the experimental data appear to be largely asymmetric for the arithmetic diameter profile D10 on the the left of fig. 10, which raises questions regarding the accuracy of experimental particle size measurement. Indeed it was found difficult to acquire particle data even close to the atomizer nozzle for the cold operating point as the spray is dense and particle evaporation is negligible. The particle size distribution should be further investigated but, considering the measurement uncertainties and the good prediction of the velocity, the simulation is estimated in sufficiently good agreement with experiments to proceed with ignition simulations. Fig. 10: Comparison of the arithmetic mean (D10) and Sauter mean (D32) diameters of the dispersed phase at 6 mm from the combustion chamber inlet plane. The experimental data is represented with blue plain circles [21], simulation results of the CEDRE platform (present work) appear in plain black lines. Numerical diameter profiles are not provided in the work of Eyssartier [START_REF] Eyssartier | Study and modelisation of stationnary and transient two-phase flow combustion[END_REF] as monodisperse Euler-Euler simulations with a particle diameter d p = 60µm were performed.

Simulation of an ignition phase

After the validation of the non-reactive flow, the LES of an ignition sequence of the Mercato configuration is now simulated using the numerical parameters presented in sections 4.1 and 4.3.2 in combination with the ignition model presented in section 2.

Turbulent combustion model

In the present simulations, the chemical kinetics of kerosene rely on the BFER scheme of Franzelli et al. [START_REF] Franzelli | A two-step chemical scheme for Large-Eddy Simulation of kerosene-air flames[END_REF]. The dynamically Thickened Flame Model (TF) designed for partially premixed turbulent combustion by Colin et al. [START_REF] Colin | A thickened flame model for large eddy simulations of turbulent premixed combustion[END_REF] is used. The sub-grid scale interactions between turbulence and the flame are modeled by the efficiency function proposed by Charlette et al. [START_REF] Charlette | A power-law wrinkling model for LES of premixed turbulent combustion. Part I -non-dynamic formulation and initial tests[END_REF]. In addition, a sensor is used to detect the position of the flame front in order to limit the thickening procedure to the vicinity of the flame.

When applying the thickened flame model to dispersed two phase flows, the source terms of the carrier phase resulting from the dispersed phase need to be divided by the thickening factor [START_REF] Boileau | Simulation aux grandes échelles de l'allumage diphasique des foyers aéronautiques[END_REF]. This consistency correction ensures that propagation speeds of laminar spray flames remain independent of the thickening factor. This correction was validated using numerical simulations of one-dimensional saturated two-phase kerosene spray flames, see fig. 11. Note that the heat release was divided by its maximum value in each simulation. Both heat release rate and droplet diameter profiles are almost unaltered by the uniform thickening procedure compared to the reference case. On the other hand, the dynamical thickening procedure appears to slightly spread the area of significant heat release, but the integral of the latter remains almost identical. The impact of the dynamical thickening procedure on the droplet evaporation rate seems even more moderate.

Numerical setup

The numerical grid was modified with respect to the previous non-reactive simulations as significant refinement was applied to the the area of energy deposition. The characteristic local mesh size was approximately set to the laminar flame thickness to guarantee that negligible thickening was applied in the first instants of the development of the flame kernel.

In order to partly reproduce the variability of ignition, the energy deposition procedure described in section 2.1 was applied at four different instants of the two-phase flow simulations. The duration of energy deposition was deduced from experimental data and fixed to 50 µs. The amount of energy supplied to the flow was 300 mJ, which is larger than what would be expected from common estimates for the current set-up as the total electric energy supplied to the spark plug amounts to 400 mJ and energy losses are estimated at about 70%, leading to an energy transfer to the flow of approximately 120 mJ. At the currently investigated spark plug location, i.e. z = 56 mm from the chamber inlet plane, ignition was experimentally never observed before the 9 th spark, which at a frequency of 6 Hz means a second and a half of fuel carburation. Therefore, the non-reactive dispersed two-phase flow simulation should have been run for the same duration in order to reproduce the experimental fuel filling of the chamber. This was however not possible for reasons of computational expense, and it was compensated by an increase in energy forcing ignition with 300 mJ after 150 ms from the start of fuel injection. At this time, the chamber seems sufficiently filled with liquid droplets to ignite. Further fueling is only expected to increase liquid volume fractions in regions where accumulation occurs, i.e. in the lateral and central recirculation zones and the chamber walls via droplet splashing. Note that although predicted by the current droplet wall interaction model, the formation of a liquid film on the chamber walls is not explictly taken into account in the present numerical simulations. Furthermore, it was verified that the liquid volume fractions in the vicinity of the spark location fluctuate around an apparently steady mean value. The increase in the minimum ignition energy when the carburation time is reduced was observed by Neophytou et al. [START_REF] Neophytou | Large-eddy simulation of ignition and flame propagation in a trisector combustor[END_REF].

Results

Among the four ignition simulations, three led to failed ignition (extinction of the kernel while it was propagating within the chamber) and one led to a successful ignition of the Mercato test rig. In the failed ignition cases, a flame kernel was created and convected by the flow but was not strong enough to grow and eventually vanished. Snapshots of the successful ignition sequence at different instants after sparking are shown on fig. 12. It appears that a flame kernel forms and is convected by the flow toward a region with an equivalence ratio close to stoichiometry at about t = 4.5 ms after energy deposition. From this point on, a propagating flame starts and develops in the chamber, the propagation becoming much faster after t = 7.5 ms. Note that in the failed ignition cases, the flame kernel remained in areas of low total equivalence ratio which tends to indicate that the equivalence ratio at the flame kernel location is a controlling parameter of ignition. In order to confirm this assumption with quantitative indicators, a flame consumption speed was evaluated using the following expression:

S c = - 1 ρ 0 Y 0 O2 A f E ωO2 F dV (9)
where ρ 0 is the fresh gas density, Y 0 O2 the oxygen mass fraction in the fresh gas and A f the flame area evaluated with the surface of an iso-contour T = 1500K. The integrated O 2 source term is the one of the Thickened Flame model, introducing the efficiency function E and the thickening factor F [START_REF] Colin | A thickened flame model for large eddy simulations of turbulent premixed combustion[END_REF]. The temporal evolution of this flame speed is displayed on fig. 13. Besides two peak values corresponding to a vanishingly small flame area, the flame speed is contained within plausible ranges for a moderately turbulent twophase kerosene flame, i.e. between 0.2 and 0.5 m s -1 . The expansion of the flame kernel is visible with the increase of the area enclosing T = 1500K iso-contour. Fig. 13: Evolution of the flame consumption speed defined according to eq. 9 and the flame area identified by a temperature isosurface T = 1500K over time.

A qualitative comparison with experiment based on flame visualisations at identical instants is provided in fig. 14. It is reminded that both ignition sequences occur with significantly different equivalence ratios in the chamber as ignition occurs at about 150 ms of fueling in the simulations compared to 1.5 s in the experiment. This important discrepancy clearly hinders detailed comparisons. Despite this fact, the previously identified ignition phases also seem present in the experiment as the reactive kernel is merely convected without growing until at least t = 7 ms. However, the flame propagation appears significantly faster in the simulations: it is already completed after 27 ms whereas it lasts 39 ms in the experiment. Nevertheless, it also seems that the propagation phase is dominated by flow convection in both cases since there is approximately an order of magnitude between the observed propagation speed of the flame kernel and the turbulent flame speed. This was already observed in Boileau [START_REF] Boileau | LES of an ignition sequence in a gas turbine engine[END_REF] and Philip [START_REF] Philip | Large Eddy Simulations of the ignition sequence of an annular multiple-injector combustor[END_REF] and confirms that hot gas expansion plays a key role in the full ignition of a burner. It may therefore be argued that some major features of ignition are captured in the present simulation.

Conclusions

The current work presented simulation results of a successful ignition sequence for the dispersed two-phase flow within the experimental test rig Mercato. Performing such ignition simulations required the implementation of simplified models describing the energy deposition process and thermophysical properties at high temperatures. These ingredients were partly validated with the simulation of a laminar freely propagating premixed flame kernel, for which present numerical results were in very good agreement with both experimental data and previous numerical results. Then, numerical simulations of the non-reacting purely gaseous and two-phase flow in the Mercato test rig were realized for validation purposes. Results for the non-reacting gaseous flow were found in excellent agreement with experiment for both mean and root mean square velocity profiles. For the two-phase flow simulations, a simplified injection model was used and proved to yield good results for velocity profiles. However, more significant differences were observed for the size distribution profiles while global levels were in reasonable agreement. Finally, ignition simulations of the Mercato configuration were performed. Analysis of the flame behavior for a successful ignition event indicated that the flame propagation could be divided in two main phases: a convection of the flame kernel into regions of low velocity and equivalence ratios close to unity where it can grow, followed by a rapid expansion within the chamber due to convection by the flow field. These findings seem confirmed by qualitative comparisons with experiment. The flame propagation is however too fast in the simulation, which may be an artefact of the different chamber fueling between the simulation and the experiment. Overall, the results obtained in the present work show that the simulation is able to reproduce qualitatively ignition sequences. Quantitative predictions require further investigation on both experimental and numerical sides. In particular, the characteristics of the two-phase mixture at the time of ignition in the experiments would be very helpful. Work is currently ongoing on the various simplified models used here, namely the spray injection and chemical kinetic scheme.

Fig. 1 :

 1 Fig. 1: Evolution of the heat capacities of N 2 over temperature according to CEDRE reference data and the data of McBride et al. [25]. Curves: numerical results from CEDRE simulations. Symbols: polynomial curves.

Fig. 2 :

 2 Fig.2: Evolution of the flame kernel radius over time for the ignition of a laminar propane-air mixture : experimental data from Renou and Boukhalfa[START_REF] Renou | An experimental study of freely propagating premixed flames at various Lewis numbers[END_REF], AVBP simulation from Enaux[START_REF] Enaux | Simulation aux Grandes Echelles d'un moteur à allumage commandéévaluations des variabilités cycliques[END_REF] and CEDRE simulation from the present work with two values of deposited energy.

Fig. 4 :Fig. 5 :

 45 Fig. 4: View of the experimental setup for the Mercato test rig

Fig. 7 :

 7 Fig. 7: Comparison of non-dimensional mean (top) and root mean square (bottom) air axial velocity profiles at several distances from the combustion chamber inlet plane. The velocities were made nondimensional by dividing them by the inlet velocity (92.45 m/s). The experimental data is represented with blue plain circles [21], simulation results of the AVBP solver [13] and CEDRE (present work) respectively appear in dashed red and plain black lines.

Fig. 8 :

 8 Fig. 8: Comparison of non-dimensional mean (top) and root mean square (bottom) tangential air velocity profiles at several distances from the combustion chamber inlet plane. The velocities were made nondimensional by dividing them by the inlet velocity (92.45 m/s). The experimental data is represented with blue plain circles [21], simulation results of the AVBP solver [13] and CEDRE (present work) respectively appear in dashed red and plain black lines.

Fig. 9 :

 9 Fig. 9: Comparison of the non-dimensional mass averaged mean axial and radial velocity profiles of the dispersed phase at 6 mm from the combustion chamber inlet plane. The velocities were made nondimensional by dividing them by the inlet velocity (92.45 m/s). The experimental data is represented with blue plain circles [21], simulation results of the AVBP solver [13] and CEDRE (present work) respectively appear in dashed red and plain black lines.

Fig. 11 :

 11 Fig. 11: Non-dimensional heat release rate and squared droplet diameter profile across a one-dimensional kerosene spray flame. Black line: reference results without thickening, i.e. F = 1. Red dashed line: uniform thickening applied to the entire computational domain (F = 5). Blue dotted line: dynamical thickening procedure based on a sensor limiting the thickening to the vicinity of the flame front (F max = 5)

Fig. 12 :

 12 Fig. 12: Numerical snapshots of an iso-contour of temperature 1500 K colored by axial velocity on a field of total equivalence ratio at different instants after energy deposition. The spark location is schematized on the top left figure (t=1 ms).

Fig. 14 :

 14 Fig. 14: Experimental tomography and simulated temperature iso-contour (T = 1500 K) snapshots at five different instants. The instants are indicated on the experimental snapshots and are identical for the simulation. Although the choice of the isosurface is arbitrary, it does not significantly affect the qualitative visual comparison.

Table 2 :

 2 Parameters of the main operating points of the Mercato test rig. The so called reference point is simulated in the present work. 4 Validation of the non-reacting flow in the Mercato configuration 4.1 Numerical setupThe numerical simulations rely on the unsteady compressible Navier-Stokes finite-volume solver CHARME of the CEDRE platform. A Large Eddy Simu-lation (LES) approach with the Smagorinsky sub-grid scale turbulence model [39] is employed. The main numerical parameters are summarized in table 3.

  [START_REF] Charlette | A power-law wrinkling model for LES of premixed turbulent combustion. Part I -non-dynamic formulation and initial tests[END_REF]. The resulting grid is composed of approximately 3,758,000 tetrahedra, which coincides with the number of degrees of freedom since CHARME relies on a cell-centred formalism. Non reflecting boundary conditions are imposed at the in-and outlet while adiabatic laws of the wall are used to model all solid boundaries. The results of the present work were averaged over approximately 450 ms, representing approximately 8 flow through times of the combustion chamber.

	Number of tetrahedra	3 758 000
	Time integration	RKI2 (2 nd order)
	Time step	∆t = 10 -6 s
	Space discretisation scheme	Roe (2 nd order)
	Turbulence model	Smagorinsky Cs = 0.1
	In-and outlet	Non-reflecting
	Solid boundaries	Law of the wall, isothermal

Table 3 :

 3 Main numerical parameters for the non-reacting simulations of the gaseous and dispersed two-phase flows

	4.2 Gaseous flow
	A comparison of simulation results with experimental data is shown on figs. 7
	and 8

Acknowledgements The financial support of the Direction Générale de l'Armement (DGA), the French Government Defense procurement and technology agency, is gratefully acknowledged. The authors would like to warmly thank Olivier Rouzaud, Lionel Matuszewski and Nicolas Bertier for useful discussions.

Compliance with ethical standards Loïc Hervo's PhD thesis was partially funded by the Direction Générale de l'Armement (DGA), the French Government Defense procurement and technology agency.

Conflict of Interest

The authors declare that they have no conflict of interest.