Nathan X Roth
email: nathaniel.x.roth@nasa.gov

Boncho Erika L Gibb

P Bonev

Michael A Disanti

Neil Dello Russo

Adam J Mckay

Ronald J Vervack

Hideyo Kawakita

Mohammad Saki

Nicolas Biver

/ 21p

Nathan X Giacobini-Zinner

Roth

Boncho P Bonev

Dominique Bockelée-Morvan

Lori Feaga

Nicolas Fougere

Anita Cochran

Michael Combi

Yinsi Shou

Probing the Evolutionary History of Comets: An Investigation of the Hypervolatiles CO, CH4, and C2H6 in the Jupiter-family Comet 21P/Giacobini-Zinner

published or not. The documents may come

L'archive ouverte

INTRODUCTION

Comets are among the most primitive remnants of the solar system's formation. They accreted in the early stages of the protosolar nebula and have been stored for the last ~4.5 Gyr in the cold outer solar system in the scattered Kuiper disk or in the Oort cloud dynamical reservoirs.

Because comets lack a known mechanism for efficient internal self-heating owing to their small sizes, the present-day volatile composition of their nuclei likely reflects to a large degree the composition and conditions where (and when) they formed. Thus, measuring the volatile composition of comets offers an opportunity to place observational constraints on the history of the early solar system by measuring the abundances of trace species in their nuclei (Mumma andCharnley 2011, Bockelée-Morvan et al. 2004).

The volatile inventory of comets can be inferred by studying their coma composition during passages into the inner solar system (heliocentric distance R h  3 au). Such studies at multiple wavelengths have revealed extensive chemical diversity among the comet population (e.g., A'Hearn et al. 1995, Crovisier et al. 2009, A'Hearn et al. 2012, Cochran et al. 2012, Ootsubo et al. 2012, Dello Russo et al. 2016). In particular, high-resolution near-infrared spectroscopy provides a valuable tool for sampling the composition of primary volatiles (i.e., ices subliming directly from the nucleus) of comets via analysis of fluorescence emission in cometary comae. To date, over 30 comets have been characterized in this manner, with the hypothesis that the primary volatile composition of the coma can be used to infer the composition of the nucleus, and can therefore be tied to nascent solar system conditions.

However, the results of recent rendezvous missions to comets, such as the Rosetta mission to comet 67P/Churyumov-Gerasimenko, have raised significant questions regarding the nature of comets. These include questions such as: How did comets form? How are comet ices put together? How do comets change with time? To what degree do comets retain cosmogonic signatures in their nuclei? How does coma composition vary throughout a perihelion passage?

(see A'Hearn 2017 for a discussion of these questions). Whereas Oort cloud comets (OCCs) can generally be observed only during a single apparition, short-period (or ecliptic) comets, primarily

Jupiter-family comets (JFCs), offer the opportunity to investigate potential evolutionary effects on volatile composition due to frequent and repeated perihelion passages, as well as a search for variability in coma composition on a variety of timescales (e.g., Knight & Schleicher 2013, Combi et al. 2019) depending on observational opportunities. Most processes that may alter the properties of the nucleus are expected to affect a thin (at most a few meters deep) layer near the surface, which is likely lost over the course of a typical perihelion passage (see Stern et al. 2003).

Nonetheless, an ecliptic comet that experiences many perihelion passages, particularly at small R h , may (potentially) experience considerable processing compared to an OCC. Indeed, measured JFCs are on average depleted in certain primary volatiles, such as C 2 H 2 and C 2 H 6 , relative to OCCs (Dello Russo et al. 2016). Understanding these potential evolutionary (or, natally inherent) effects, including observed differences between JFCs and OCCs, is essential for placing the results of coma composition studies into a meaningful context in the framework of solar system formation.

As the most volatile molecules systematically observed in comets, the "hypervolatiles" CO, CH 4 , and C 2 H 6 may be the most sensitive to both primordial conditions as well as thermal processing in comets (Dello Russo et al. 2016), and characterizing each molecule in an individual ecliptic comet may provide unique insights into its evolutionary history. Being symmetric hydrocarbons, CH 4 and C 2 H 6 can only be sampled in the near-infrared due to their lack of a dipole moment and thus lack of allowed rotational modes. Of these molecules, detection of C 2 H 6 has been reported in ten ecliptic comets to date. Whereas CO is readily detectable at radio wavelengths and has been measured in several ecliptic comets (Crovisier et al. 2009), near-infrared measurements of CO and CH 4 in such comets are more elusive. The reason for this is that measuring near-infrared transitions of CO and CH 4 requires sufficiently large geocentric velocities ( dot ) to Doppler shift cometary emission lines away from their corresponding highly opaque telluric counterparts and into regions of adequate atmospheric transmittance. However, the overall lower gas production rates and hence fainter nature of JFCs (compared with OCCs) means that most observations traditionally take place near closest approach to Earth, coinciding with insufficient  dot to measure their emission lines from CO and (especially) CH 4 . This results in a paucity of complete hypervolatile inventories for JFCs (e.g., see Dello Russo et al. 2016).

Fortunately, the increased sensitivity and long on-source integration times afforded by the high-resolution iSHELL spectrograph, which became available for use at the NASA Infrared Telescope Facility (IRTF) in late 2016, together with unusually favorable apparitions for several short-period comets have to date enabled sensitive searches for these molecules in short-period comets 2P/Encke [START_REF] Roth | [END_REF]) and JFC 45P/Honda-Mrkos-Pajdušáková (DiSanti et al.

2017

). The highly favorable 2018 perihelion passage of JFC 21P/Giacobini-Zinner (hereafter G-Z) featured sufficiently high geocentric velocity simultaneously with small geocentric distance, and afforded the opportunity to characterize its hypervolatile content. G-Z is the prototype for the eponymous "GZ-type" of carbon-chain depleted comets, depleted in both C 2 and NH 2 with respect to H 2 O and accounting for ~6% of comets measured (A'Hearn et al. 1995, Fink 2009).

Coupled with published near-infrared observations of G-Z during the 1998 and 2005 apparitions (Weaver et al. 1999, Mumma et al. 2000, DiSanti et al. 2013), our measurements also enabled searches for coma compositional variability of hypervolatiles on timescales of days, both pre-vs.

post-perihelion, and across multiple perihelion passages. In Section 2, we discuss our observations and data analysis. In Sections 3 and 4, we present our results. In Section 5, we compare our results to those from previous perihelion passages. In Section 6, we examine G-Z's place in the context of other comets characterized to date.

OBSERVATIONS AND DATA REDUCTION

During its 2018 apparition, G-Z both reached perihelion (1.01 au) and was closest to Earth (0.39 au) on September 10. On UT 2018 July 25, 28, 29, and 31, September 7 and 11, and October 10, we observed G-Z with the high-resolution (λ/Δλ ~ 40,000), near-infrared, immersion-grating echelle spectrograph iSHELL [START_REF] Rayner | Proc. SPIE[END_REF][START_REF] Rayner | Proc. SPIE[END_REF] at the 3-m NASA IRTF to characterize its hypervolatile composition. We utilized three iSHELL settings (Lcustom, Lp1, and M2) so as to fully sample a suite of molecular abundances. We oriented the slit along the projected Sun-comet line on all dates (see Table 1).

On our July and October dates, observations were performed with a 6-pixel (0.75) wide slit, using our standard ABBA nod pattern, with A and B beams symmetrically placed about the midpoint along the 15 long slit and separated by half its length. A malfunction of the iSHELL dekker precluded the use of the 15 long slit in September and necessitated off-chip nodding using the 0.75 wide by 5 long slit. We placed G-Z in the center of the slit for A frames and nodded 20 perpendicular to the slit for B (sky) frames. Combining spectra of the nodded beams as A-B-B+A cancelled emissions from thermal background, instrumental biases, and "sky" emission (lines and continuum) to second order in air mass. Flux calibration was performed using appropriately placed bright IR flux standard stars on each date using a wide (4.0) slit. On October 10, observing time lost owing to a telescope pointing error precluded the acquisition of flux calibration sets; therefore, we adopted a flux calibration factor ( W/m 2 /cm -1 /counts/s), based on that measured on other dates. Although this could affect absolute production rates, derived mixing ratios should be unaffected, as our targeted molecules were observed simultaneously or contemporaneously with water or OH prompt emission (OH*, a proxy for water production, see Bonev et al. 2006). The observing log is shown in Table 1.

[Insert Table 1 Here]

Our data reduction procedures have been rigorously tested and are described extensively in the refereed literature (Bonev 2005, DiSanti et al. 2006, Villanueva et al. 2009, Radeva et al. 2010, DiSanti et al. 2014), and their application to unique aspects of iSHELL spectra is detailed in §3.2 of DiSanti et al. (2017). Here we will only discuss aspects of the reduction of G-Z frames which differed from those previously reported.

Each echelle order within an iSHELL setting was processed individually as previously described, such that each row corresponded to a unique position along the slit, and each column to a unique wavelength. We found that spatially resampling using a third-order polynomial more completely removed the curvature in the spatial dimension from iSHELL frames, and so employed this in place of previously used second-order polynomial resampling (DiSanti et al. 2017). Spectra were extracted from the processed frames by summing signal over 15 rows, seven rows to each side of the nucleus, defined as the peak of dust emission in a given spectral order.

For our September observations (using the shorter 5 slit), we found that the iSHELL flat lamp provided illumination of the chip that was not consistent with that of the sky. This introduced a curvature effect into our spectra, which we corrected by fitting and then subtracting a polynomial baseline (Fig. 1). This may have affected flux calibration and therefore the calculated absolute molecular production rates (Q's) reported for our September dates. However, emissions from all molecules within each individual iSHELL setting are sampled simultaneously, and therefore our derived mixing ratios should be unaffected. The general consistency of production rates and derived mixing ratios on both September dates suggests that any uncertainties introduced by this illumination offset were likely minimal.

[Insert Figure 1]

We determined contributions from continuum and gaseous emissions in our comet spectra as previously described (e.g., DiSanti et al. 2016), and illustrate the procedure in Figure 2. We convolved the fully resolved transmittance function to the resolving power of the data (~ 4.0 x 10 4 ) and scaled it to the level of the comet continuum. We then subtracted the modeled continuum to isolate cometary emission lines and compared synthetic models of fluorescence emission for each targeted species to observed line intensities.

[Insert Figure 2] Nucleocentric (or "nucleus-centered") production rates (Q NC ) were determined using our well-documented formalism (Dello Russo et al. 1998, DiSanti et al. 2001, Bonev 2005, Villanueva et multiplied by an appropriate growth factor (GF), determined using our Q-curve methodology (e.g., Dello Russo et al. 1998, DiSanti et al. 2001, Bonev 2005, Gibb et al. 2012), to establish total (or global) production rates (Q). This GF corrects for atmospheric seeing, which suppresses signal along lines of sight passing close to the nucleus due to the use of a narrow slit, as well as potential drift of the comet during an exposure sequence. Global production rates for all detected molecules are listed in Table 2. GFs were determined for both gas and dust when the signal-tonoise ratio (S/N) was sufficiently high (i.e., only for H 2 O, CO, and C 2 H 6 ). For September dates, the short 5 slit precluded the use of Q-curves to calculate GFs. Therefore, we assumed a GF of 1.8, a value consistent with that obtained from our other dates (Table 2).

RESULTS

Spatial Profiles

For July dates, we were able to extract spatial profiles for H 2 O, CO, and C 2 H 6 emissions in G-Z (Fig. 3). Within uncertainty, this suggests that emission for all three species closely tracked that of the co-measured dust; therefore, on dates for which molecular GFs could not be well-constrained, we adopted that of dust co-measured within each setting when calculating production rates (Table 2).

[Insert Figure 3 Here]

Mixing Ratios of Volatile Species

Molecular Fluorescence Analysis

Synthetic models of fluorescence emission for each targeted species were compared to observed line intensities, after correcting each modeled line intensity (g-factor) for the monochromatic atmospheric transmittance at its Doppler-shifted wavelength (according to the geocentric velocity of the comet at the time of the observations). The g-factors used in synthetic fluorescence emission models in this study were generated with quantum mechanical models developed for CO (Paganini et al. 2013), CH 4 (Gibb et al. 2003), C 2 H 6 (Villanueva et al. 2011b), and H 2 O (Villanueva et al. 2012). A Levenberg-Marquardt nonlinear minimization technique (Villanueva et al. 2008) was used to fit fluorescence emission from all species simultaneously in each echelle order, allowing for high-precision results, even in spectrally crowded regions containing many spectral lines within a single instrumental resolution element. Production rates for each sampled species were determined from the appropriate fluorescence model at the rotational temperature of each molecule ( §3.2.2).

Determination of Rotational Temperature

Rotational temperatures (T rot ) were determined using correlation and excitation analyses as described in Bonev (2005, pp. 53-65), Bonev et al. (2008), DiSanti et al. (2006), and Villanueva et al. (2008). In general, well-constrained rotational temperatures can be determined for individual species having intrinsically bright lines and for which a sufficiently broad range of excitation energies is sampled. Utilizing the large spectral grasp of iSHELL, in the case of H 2 O we were able to sample dozens of strong lines simultaneously.

We found consistent rotational temperatures for multiple molecules on all dates (including H 2 O). The T rot for H 2 O was well-constrained on September 7 (being 75 ± 3 K) and was consistent (within 1σ uncertainty) with that for C 2 H 6 on September 11 (66 -9 +12 K). Rotational temperatures for our July dates were also in formal agreement, being 64 -11 +15 K for CO on July 28 and 48 -13 +19 K for H 2 O on July 29. We were unable to derive well-constrained rotational temperatures for any molecules on October 10. We calculated production rates and mixing ratios at T rot = 48 K and 64 K for July dates and varied T rot as a parameter for the October date, calculating production rates and mixing ratios for each molecule at representative values T rot = 50 K, 60 K, and 70 K. In general, mixing ratios for a given species derived at each temperature are consistent with one another within 1σ uncertainty (Table 2).

[Insert Table 2 Here]

Secure Detections of Hypervolatiles

Our detections of CO, CH 4 , and C 2 H 6 in G-Z are particularly notable for two reasons: (1)

They address the paucity of measurements of CO and CH 4 in ecliptic comets in general, and (2)

The measurement of these hypervolatiles in an individual ecliptic comet across multiple perihelion passages, and on multiple dates during its 2018 apparition. Of all primary volatiles systematically measured in comets, these three molecules are most sensitive to thermal processing, but as noted earlier, CO and CH 4 are also among the most difficult to sample from the ground due to lack of sensitivity and/or adequate geocentric velocity. G-Z's excellent geocentric velocity (| dot | > 13 km s -1 pre-perihelion, | dot | > 11 km s -1 post-perihelion) allowed firm detections of all three species. CO and CH 4 have been measured in fewer than ten ecliptic comets (with most detections being below the 5σ level), making our measurements in G-Z a critical component in establishing statistics for these species in ecliptic comets, and determining the importance of natal versus evolutionary effects on present cometary volatile composition.

Coma Hypervolatile Composition Throughout the 2018 Perihelion Passage of G-Z

The 2018 apparition of G-Z provided an opportunity to conduct the first comprehensive comparison of hypervolatile abundances for a comet through three perihelion passages and also on multiple dates within a given perihelion passage, thereby allowing us to address pressing questions in cometary science. These include testing possible evolutionary effects on coma volatile composition, as well as searching for coma compositional variability on multiple timescales, including day-to-day, pre-vs. post-perihelion, and across perihelion passages. We discuss each of these topics in turn, and place G-Z in the context of other comets observed to date. Unless otherwise noted, all dates refer to the 2018 apparition.

CO

We found clear, simultaneously measured detections of CO and H 2 O on multiple dates in G-Z (Figs. 2,4(A),4(D)) pre-as well as post-perihelion. The mixing ratio CO/H 2 O was consistent on all pre-perihelion July dates within 1σ uncertainty (Table 2; see also Fig. 5) with a weighted average abundance of 1.72 ± 0.12% for T rot = 64 K. This was somewhat lower post-perihelion in October, 1.26 ± 0.23% for T rot = 60 K, suggesting that CO/H 2 O in G-Z may display pre-vs. postperihelion asymmetry. However, given the uncertainty in T rot in October, it is important to note that the range of possible October CO mixing ratios is in formal agreement with those from July (Table 2). These mixing ratios are depleted with respect to the mean for all comets measured to date at near-infrared wavelengths (5. and derived a meaningful constraint on its mixing ratio on October 10 (see Fig. 4(B) and Table 2). This suggests that CH 4 may have been variable from day-to-day in G-Z. However, there are important caveats for our CH 4 study.

We were unable to derive a well-constrained rotational temperature for CH 4 owing to the small spread in excitation energies in the sampled lines (see §3.2.2, Fig. 4B). However, we found that calculated CH 4 production rates and mixing ratios showed a sensitive dependence on assumed T rot (

C 2 H 6

Of all the hypervolatiles, C 2 H 6 is the most routinely sampled in comets owing to its intrinsically strong near-infrared transitions and the availability of multiple emissions in regions of favorable telluric transmittance independent of  dot . This enabled us to measure C 2 H 6 mixing ratios on multiple dates during G-Z's 2018 apparition, including pre-perihelion, near perihelion, and post-perihelion. Similar to CO, we found C 2 H 6 mixing ratios were consistent (within uncertainties) pre-perihelion (weighted average 0.24 ± 0.03% for T rot = 64 K) and additionally near perihelion (weighted average 0.29 ± 0.02%). However, C 2 H 6 was lower post-perihelion with mixing ratio 0.14 ± 0.03% (assuming T rot = 60 K). Compared to ecliptic comets measured to date, pre-perihelion and near perihelion G-Z was consistent with the mean mixing ratio measured for C 2 H 6 (0.34 ± 0.07%), but was depleted post-perihelion (and was severely depleted compared to the mean for all comets measured (0.55 ± 0.08%), including ecliptic comets). In the same manner as CH 4 , the use of Q(H 2 O) from H 2 O in the M2 setting rather than from OH* in the Lp1 setting to calculate C 2 H 6 mixing ratios in July and October may have contributed to its suggested variability. Similarly, the use of assumed GFs for C 2 H 6 on some of our dates may have introduced additional uncertainty into the mixing ratio C 2 H 6 /H 2 O. However, we note that Q(H 2 O) was dramatically lower in October compared with both July and September, being closer to 2 x 10 28 than to 3 x 10 28 molecules s -1 (see Table 2), consistent with asymmetry in water production with respect to perihelion found by A' Hearn et al. (1995). This in turn could indicate distinct regions of the nucleus dominating the activity in G-Z at these times, and that its chemical composition (at least in terms of CH 4 /H 2 O and C 2 H 6 /H 2 O) may also be different on October 10 compared with our pre-/near-perihelion dates. 3 gives a similar comparison numerically. Figure 5 and Table 3 suggest that each hypervolatile may display at least some degree of variability, whether across perihelion passages or during a particular apparition. We discuss each species in turn.

Comparison with Previous

[Insert Figures 5A-B Here]

[Insert Table 3 Here]

CO

In the case of CO, our pre-and post-perihelion mixing ratios are consistent with (yet lower than) the upper limit from the 1998 apparition (using CSHELL) found by Weaver et al. (1999), but considerably lower than the mixing ratio reported by Mumma et al. (2000) In the context of preserving natal solar system signatures in the nucleus ices of JFCs, it is important to note the stark contrast of our CO measurements along with those reported by than on earlier 2018 dates, yet consistent with the 2005 pre-perihelion value, suggesting possible short-term (i.e., "diurnal", perhaps associated with nucleus rotation, or seasonal, such as that seen by Rosetta at 67P/Churyumov-Gerasimenko; see §5.2) variability in its C 2 H 6 abundance ratio when compared with our pre-perihelion results. It is important to note that the possible variability in C 2 H 6 implied by our G-Z measurements (as well as those from previous perihelion passages) is small compared to the overall spread of C 2 H 6 abundances in all comets measured (Fig. 5).

Discussion of Possible Variability of Coma Hypervolatile Abundances in G-Z

Combined with previous work, our results suggest that coma hypervolatile abundances in G-Z may be variable. At 67P/Churyumov-Gerasimenko, the Rosetta mission found that nucleus shape and the location of active areas, combined with seasonal and rotational illumination effects, resulted in coma compositional variability on a variety of timescales. To further examine the nature of the suggested variability in G-Z, we examined the evolution of molecular production for each species reported here during the 2018 apparition.

Figure 6 shows the production rate of each species relative to perihelion (on September 10). Our measurements for all four species (H 2 O, CO, CH 4 , C 2 H 6 ) are consistent with A'Hearn et al.

(1995), who found that G-Z was more active pre-perihelion than post-perihelion. However, our results indicate that the relative asymmetry in molecular production is more pronounced for the trace species than for H 2 O, which is reflected in their generally lower post-perihelion compared to pre-perihelion mixing ratios (Fig. 5, Table 3).

[Insert Figure 6 Here]

In order to test whether the possible variability indicated by our results is owing to the activity of H 2 O vs. that of the trace species in G-Z, we compared the ratios CO/C 2 H 6 from the 2018 perihelion passage. We found that CO/C 2 H 6 was consistent within uncertainty pre-vs. post- It is important to note that the variability suggested by the measurements in Table 3 cannot explain the variation in each molecule among comets revealed in measurements to date (Figure 5). If the range of volatile abundances observed among all comets can be reproduced by time-resolved observations of one comet, we could seriously question the extent to which chemical diversity among the population is cosmogonic. Alternatively, a comet such as G-Z, in which measurements over three perihelion passages suggest (with carefully explored caveats) that the abundances of CO/H 2 O and C 2 H 6 /H 2 O vary on scales much smaller than the comet-tocomet range, may serve as evidence that we are indeed sampling cosmogonic signatures in our present-day measurements of parent volatiles in short-period comets. Further unraveling the complex relationship between nascent solar system conditions and evolutionary processes in comets clearly requires increasing the sample size of serial measurements in short-period comets, particularly observations targeting hypervolatiles.

Comparison to Comets as Measured at Near-Infrared Wavelengths

Comprehensive hypervolatile abundances have been securely measured in 19 OCCs to date, yet in only four ecliptic comets, including G-Z. This highlights that statistics for these species in ecliptic comets (particularly CO and CH 4 ) are far from being firmly established.

Figure 7 shows relative hypervolatile abundances reported in all comets to date, including G-Z and measurements taken by Rosetta at 67P/Churyumov-Gerasimenko using ROSINA The relative isolation of G-Z in Figure 7 further highlights the spread in hypervolatile abundances among ecliptic comets, reflecting its unique combination of CH 4 being consistent with the mean among OCCs versus the moderately depleted values for CO and C 2 H 6 . This underscores that much work remains in firmly characterizing the ranges of hypervolatile abundances in ecliptic comets and understanding their implications for placing such measurements into a meaningful context.

Summary of Results

We characterized the hypervolatile composition of the prototypical "GZ-type" comet 21P/Giacobini-Zinner with the powerful, recently commissioned iSHELL spectrograph at the NASA-IRTF on four pre-perihelion dates, two dates near perihelion, and one post-perihelion date. Combined with previous work, our results suggest that coma abundances of all three hypervolatiles (CO, CH 4 , and C 2 H 6 ) may be variable on several timescales, including day-to-day, pre-vs post-perihelion, and even across perihelion passages. However, as noted in §5, there are important caveats to our study, and additional serial measurements of G-Z are needed to confirm possible variability in its coma hypervolatile content. In any case, our results suggest that:

1. Mixing ratios of CO were consistent (within uncertainty) day-to-day preperihelion, but were slightly lower post-perihelion. Our measurements are consistent with depleted values compared to the mean among measured comets, as well as with an upper limit reported from the 1998 perihelion passage (Weaver et al. 1999).

2. Our measurements of CH 4 , the most severely underrepresented hypervolatile in studies of ecliptic comets, represent its first reported values in G-Z. CH 4 abundances were consistent with mean values among all comets measured, and may have been variable from day-to-day. However, there are important caveats to the possible variability of CH 4 in G-Z ( §4.2, 5.2).

3. We found that the mixing ratio of C 2 H 6 decreased significantly pre-vs. postperihelion, its post-perihelion value being consistent with being strongly depleted.

Our pre-perihelion C 2 H 6 mixing ratios were enriched compared to measurements during the same seasonal phase in 2005 (DiSanti et al. 2013), yet our postperihelion mixing ratio was consistent with the result from 2005.

4. If G-Z is indeed variable, the spread among our measurements, as well as between those from previous perihelion passages, is significantly smaller than the variation in each molecule among all comets measured (Fig. 5). This may be evidence that natal conditions dominate over evolutionary effects due to successive perihelion passages in setting the composition of short-period comets.

Understanding the cause(s) of the considerable spread of hypervolatile abundances among comets (both OCCs and short-period comets) seen in Figure 7 is necessary for disentangling primordial from evolutionary effects in setting the present-day (observed) abundances of hypervolatiles (and of primary volatiles in general) in comets. On the one hand, chemical models of protoplanetary disks (e.g., Drozdovskaya et al. 2016, Willacy et al. 2015) predict that comets incorporated a wide range of hypervolatile abundances from their formation region(s) in the protosolar nebula. On the other hand, the nontrivial effects of heterogeneous outgassing and seasonal illumination on coma composition, such as that seen by Rosetta at 67P (i.e., Bockelée-Morvan et al. 2016, Hässig et al. 2015, Luspay-Kuti et al. 2015, Fougere et al. 2016aand b, Feldman et al. 2018), cannot be overlooked. This emphasizes the high impact of serial observations of comets, particularly those targeting hypervolatiles in ecliptic comets, which may be most indicative of the role that primordial vs. evolutionary effects play in setting the composition of comets. Thankfully the availability of next-generation instruments such as iSHELL, capable of delivering the long on-source integration times and excellent sensitivity required for such measurements, is enabling us to better understand the interplay between nascent solar system conditions, evolutionary processing, and coma compositional variability when interpreting the results of primary volatile studies in comets. Weaver et al. 1999;orange, Mumma et al. 2000), 2005(pink, DiSanti et al. 2013), and 2018 (green, yellow, cyan;this work) c Global production rate. Uncertainties in production rate include line-by-line deviation between modeled and observed intensities and photon noise (see Bonev, 2005;Dello Russo et al. 2004;Bonev et al. 2007).

d Molecular abundance with respect to H 2 O.

e Continuum (dust) growth factor.

f A growth factor for September dates could not be derived owing to use of the 5 short slit; therefore, a GF of 1.8, consistent with growth factors derived for species pre-and post-perihelion, was assumed. e Mean values and 1σ uncertainties among measured comets taken from Dello Russo et al. (2016). The number of measurements used to calculate the mean is given in parentheses. Mean values for CO and CH 4 are given for OCCs only owing to the extreme paucity of such measurements in JFCs, whereas the mean for C 2 H 6 is given for all comets measured (JFCs and OCCs).

f Range among comets measured after Dello Russo et al. 2016, DiSanti et al. 2017, Roth et al. 2017[START_REF] Roth | [END_REF], Faggi et al. 2018, Dello Russo et al. 2019.

Figures 4 (

 4 Figures 4(A)-(E) show clear CO, H 2 O, CH 4 , C 2 H 6 , and OH* emissions in G-Z superimposed on

  Perihelion Passages5.1 Comparison with Previous Perihelion Passages of G-ZG-Z is the only comet observed at near-infrared wavelengths during three different perihelion passages: 1998(Weaver et al. 1999, Mumma et al. 2000), 2005(DiSanti et al. 2013), and 2018 (this work), and is just the second comet to have a comprehensive comparison of hypervolatile abundances across apparitions (the other being 2P/Encke; see[START_REF] Roth | [END_REF], Radeva et al. 2010). Figure5Ashows our individual 2018 measurements of hypervolatile abundances in G-Z. Figures 5B shows our mean pre-, near-, and post-perihelion hypervolatile abundances together with those for G-Z from 1998 and 2005, and (for comparison) those for all measured comets. Table

Figure 1 .

 1 Figure 1. Extracted spectra of comet G-Z showing order 155 of the iSHELL Lp filter taken with the 5"

Figure 2 .

 2 Figure 2. Extracted spectra showing clear detections of CO and H 2 O in comet G-Z superimposed on the cometary continuum on UT 2018 July 28. The gold trace overplotted on the uppermost cometary

Figure 3 .

 3 Figure 3. A. Spatial profiles of co-measured emissions in G-Z for H 2 O (black), CO (orange) and dust (red) on UT 2018 July 29. The slit was oriented along the projected Sun-comet line (position angle 223°),

Figure 4 .

 4 Figure 4. A. Extracted spectra showing detections of CO and H 2 O in comet G-Z on UT 2018 July 28, with traces and labels as described in Figure 2. B. Detections of CH 4 , C 2 H 6 , CH 3 OH, and OH* (prompt emission) on UT 2018 July 31. The zoomed subplots highlight the locations of individual (observed and

Figure 5 .

 5 Figure 5. A. Comparison of mixing ratios (abundances relative to H 2 O, expressed in %) of hypervolatiles in G-Z sampled on each date during the 2018 apparition. B. Comparison of mixing

Figure 6 .

 6 Figure 6. Evolution of molecular production in G-Z throughout the 2018 perihelion passage for H 2 O (left panel) and CO, CH 4 , and C 2 H 6 (right panel) with respect to perihelion (UT 2018 September 10). Error bars indicate measurements, whereas downward arrows indicate 3σ upper limits.

Figure 7 .

 7 Figure 7. Abundances ratios of hypervolatiles in comets characterized to date, including comets G-Z (this

  2 ± 1.3%), but are consistent with the few measurements in ecliptic comets(Dello Russo et al. 2016, DiSanti et al. 2017, Roth et al. 2018, Dello Russo et al.

	2019).
	4.2 CH 4

CH 4 bears the distinction of being the most severely undersampled hypervolatile in ecliptic comets, having been firmly measured in only six to date (Dello

Russo et al. 2016 and refs. therein, DiSanti et al. 2017[START_REF] Roth | [END_REF], Dello Russo et al. 2019)

. Utilizing the large spectral grasp of iSHELL, we detected CH 4 in G-Z at the 4σ level on July 28, at >6σ on July 31,

Table 2

 2 Regardless of which T rot we adopt, G-Z is not the first instance of a CH 4 -enriched JFC, with similar mixing ratios reported in 45P/Honda-Mrkos-Pajdušáková(DiSanti et al. 2017, Dello Russo et al. 2019). Our (3σ) upper limit for October 10 is similarly sensitive to assumed T rot (Table2), but is consistent with our July measurements assuming T rot  60 K, a reasonable assumption given the rotational temperatures measured in July and September for other molecules at similar R h .Additionally, we were unable to extract spatial profiles for CH 4 emission due to low

	). Assuming T rot = 48 K (from H 2 O on July 29), pre-perihelion CH 4 mixing
	ratios in G-Z (0.63 ± 0.17% for July 28, 1.12 ± 0.26% for July 31) are consistent with mean
	values in measured OCCs (0.88 ± 0.10%), yet are enriched compared to the few measurements in
	JFCs. Adopting a higher T

rot (e.g., 64 K from CO on July 28) increases the degree of CH 4 -enrichment.

S/N along the slit; therefore, we assumed GFs measured from other species (or co-measured dust) within a given date for CH 4 in order to calculate global production rates. It is possible that CH 4 outgassing differed day-to-day from that for H 2 O, CO, or co-measured dust and that the suggested variability may be due to our assumed GFs for CH 4 . That being said, we did not find any unusual outgassing patterns among the other molecules or dust relative to one another in G-Z, so we expect our assumed GF for CH 4 to be reasonable.

Finally, OH* was weak in G-Z for our July and October dates, and was only firmly detected near perihelion in September.

We therefore calculated mixing ratios for CH 4 using Q(H 2 O) obtained from the M2 setting on the same date. We estimate the inter-setting calibration uncertainty to be ~10%, and have incorporated this into the reported uncertainty in our mixing ratios. Use of contemporaneously (but not simultaneously) measured Q(H 2 O) for CH 4 abundances in July and October may account for some of the spread in abundances from date to date. However, the formal agreement between Q(H 2 O) obtained from OH* (in Lp1) and H 2 O (in Lcustom) on September 7 (Table 2) argues against both a systematic difference in retrieving water production rates in these two ways, and also against short-term variations in Q(H 2

  Clearly, further measurements of G-Z are necessary to clarify the possible variability of its coma CH 4 content.

O) in G-Z.

  from observations conducted approximately three weeks earlier. However,Mumma et al. did not detect H 2 O, even though the strong line near 2151 cm -1 (as we show in Figs.4A and 4D) was encompassed together with the CO R0 and R1 lines in the same CSHELL setting. Instead, their value for CO/H 2 O was inferred from the measured CO abundance relative to C 2 H 6 (detected at the ~5σ confidence level), and an adopted (1.9σ) value for Q(H 2 O) based on residual flux at the Doppler-shifted frequency of the 2151 cm -1 line. In any case, results obtained to date suggest that the abundance ratio of CO in G-Z may display variability, both during a single apparition (as is also suggested by our 2018 measurements) and across multiple apparitions.It is important to note that the 1998 measurements of G-Z with CSHELLthe small spectral grasp of which precluded measuring H 2 O simultaneously with CH 4 or C 2 H 6introduced uncertainties due to inter-setting calibration in addition to potential temporal variations in production rates. In contrast, the large spectral grasp of iSHELL enabled simultaneous measurements of all three hypervolatiles with either H 2 O or OH* during the 2018 perihelion passage of G-Z.

however, this mixing ratio suffers the same systematic uncertainty noted for CO/H 2

  <2σ) of H 2 O. Nonetheless, our measurements do not support this conclusion, and instead indicate that G-Z has a CO abundance that is more similar to the few measurements in ecliptic comets and is depleted when compared to all comets measured.

	variable on timescales of days to months in 2018. However, as previously noted, there are
	important caveats regarding its purported variability. In any case, our measurements indicate that
	G-Z is consistent with to enriched compared to the mean CH 4 abundance for all comets
	measured.
	5.1.3 C 2 H 6
	Our C 2 H 6 mixing ratios obtained pre-perihelion and near perihelion were consistent with
	that found by Mumma et al. (2000) (~0.2% relative to H 2 O;
	5.1.2 CH 4
	CH 4 has not been reported previously in G-Z -for both the 1998 and 2005 observations

Weaver et al. compared to those of Mumma et al. The CO/H 2 O mixing ratio inferred by Mumma et al. (10 ± 6%) would place G-Z as the only known CO-enriched JFC to date. If G-Z were

indeed so enriched in CO, it would have profound implications for the origins and evolutionary processing history of JFCs. However, as mentioned previously this is based on an extremely tentative "detection" (O in §5.1.1), but were significantly higher than the upper limit (<0.05-0.08%) reported by

Weaver et al. (1999)

, and also the measurement from 2005 (0.14%,

DiSanti et al. 2013)

. However, our post-perihelion measurement for C 2 H 6 on October 10 was considerably lower (by approximately a factor of 2)

  Hässig et al. (2015) found long-term variation in the coma abundances of CO and CO 2 due to seasonal illumination effects; furthermore, other species (such as CH 4 ) varied on shorter timescales, showing diurnal variations that differed from those of other volatiles, such as CO and C 2 H 6 (Luspay-Kuti et al.

	2015, Bockelée-Morvan et al. 2016, Fink et al. 2016). Similar effects may have contributed to
	the suggested coma hypervolatile variability in G-Z. Unfortunately, our ground-based
	measurements do not have sufficient spatial resolution to test this possibility.

differences in observational circumstances, techniques, and analysis must be kept in mind when comparing results from studies of different comets.

  measurements of C 2 H 6 (LeRoy et al. 2015), MIRO measurements of CO(Biver et al. 2019), and VIRTIS measurements of CH 4 (Bockelée-Morvan et al. 2016). The particularly low C 2 H 6 /CH 4 Similarly, as the hypervolatile abundances of more ecliptic comets are completely characterized, it appears that they may span a similar range of CO/CH 4 and C 2 H 6 /CH 4 as that observed among OCCs, from severely depleted (45P/ Honda-Mrkos-Pajdušáková #20 in Fig. 7; DiSanti et al. 2017) to near-mean values (2P/Encke, #10; Roth et al. 2018) to (possibly) enriched values in 67P(LeRoy et al. 2015, Bockelée-Morvan et al. 2016, Biver et al. 2019). It is important to note that comparisons between the in-situ measurements of LeRoy et al. (2015) and bulk coma measurements (e.g., this work, Bockelée-Morvan et al. 2016, Biver et al. 2019) are not straightforward, particularly given the differences in R h between each set of measurements. Additionally,

	2 H 6 /CH 4
	supports this hypothesis, and together with the observed blueish polarization, may indicate that
	simple hydrocarbons were efficiently converted into more complex organics in the ices that were
	incorporated into the nucleus of G-Z.
	[Insert Figure 7 Here]
	It is apparent from Figure 7 that hypervolatile abundances among OCCs span a large
	range of values.

ratio in G-Z (points 22-23) is supported by observations at other wavelengths.

Kiselev et al. (2000)

reported a blueish linear polarization spectrum for continuum in G-Z at optical wavelengths, and suggested this was caused by the presence of organic grains (or large-sized complex organics). This implies that G-Z is depleted in simple organics, such as C 2 H 6 , but is enriched in more complex organics, which may indicate warmer conditions were present during the formation and subsequent evolution of G-Z's constituent ices. Our low measured C

Measurements shown from 2018 are given as weighted averages for pre- perihelion and perihelion dates, and assume T rot = 64 K for July, T rot = 75 K for September, and T rot = 60 K for October.

  , as well as near-infrared measurements of each volatile in OCCs (blue) and ecliptic comets (red) measured to date, and the respective mean values for CO and CH 4 among OCCs and for C 2 H 6 among all comets (black, DelloRusso et al. 2016, DiSanti et al. 2017, Roth et al. 2017[START_REF] Roth | [END_REF], Faggi et al. 2018, Dello Russo et al. 2019). Error bars indicate measurements with 1 uncertainties, whereas downward arrows indicate 3σ upper limits.

  detected.1 Values for 67P using C 2 H 6 /H 2 O as reported in LeRoy et al. (2015) for the northern hemisphere. 2 Values for 67P using C 2 H 6 /H 2 O as reported in LeRoy et al. (2015) for the southern hemisphere. In both cases, the blue ovals trace the total possible spread in 67P's hypervolatile content.

Table 1

 1 Observing Log for 21P/Giacobini-Zinner

	UT Date	iSHELL	UT	R h (au) dR h /dt	Δ (au)	dΔ/dt	T int	Slit PA/
	(2018)	Setting			(km s -1 )		(km s -1 )	(minutes)	Length
	July 25	M2	12:02-13:58	1.20	-12.72	0.64	-13.67

Table 2

 2 Hypervolatile Composition of Comet 21P/Giacobini-Zinner

	iSHELL	Molecule	T rot	a	GF b	Q c	Q x /Q H2O	d
	Setting		(K)		(10 25 mol s -1 )	(%)
		2018 July 25, R h = 1.20 au, Δ = 0.64 au, dΔ/dt = -13.67 km s -1
	M2	H 2 O	(48)	1.82 ± 0.17 e	2692 ± 292	100
		CO	(48)	(1.82)	40.5 ± 5.1	1.51 ± 0.25
		H 2 O	(64)	(1.82)	3028 ± 306	100
		CO	(64)	(1.82)	47.4 ± 5.9	1.56 ± 0.25
		2018 July 28, R h = 1.18 au, Δ = 0.61 au, dΔ/dt = -13.5 km s -1
	Lp1	C 2 H 6	(48)	(1.91)	6.26 ± 1.28	0.23 ± 0.05
		CH 4	(48)	(1.91)	17.5 ± 4.2	0.63 ± 0.17
		C 2 H 6	(64)	(1.91)	7.02 ± 1.33	0.24 ± 0.05
		CH 4	(64)	(1.91)	26.1 ± 6.2	0.88 ± 0.24
	M2	H 2 O	(48)	1.91 ± 0.14 e	2771 ± 251	100
		CO	(48)	(1.91)	45.1 ± 3.6	1.63 ± 0.20
		H 2 O	(64)	(1.91)	2961 ± 297	100
		CO	64 -11 +15	(1.91)	50.4 ± 4.0	1.70 ± 0.22
		2018 July 29, R h = 1.17 au, Δ = 0.61 au, dΔ/dt = -13.3 km s -1
	Lcustom	H 2 O	48 -13 +19	(1.97)	2643 ± 229	100
	M2	H 2 O	(48)	(1.97)	2527 ± 345	100
		CO	(48)	1.97 ± 0.21	34.8 ± 4.7	1.38 ± 0.26

a Rotational temperature. Values in parentheses are assumed.

b Growth factor. Values in parentheses are assumed.

Table 3

 3 Hypervolatile Abundances in 21P/Giacobini-Zinner Across ApparitionsUpper limits for non-detected species are 3σ. In all cases values are expressed as % relative to H 2 O. Abundances for CO and C 2 H 6 are given as weighted averages for molecules detected on multiple dates, assuming T rot = 64 K for pre-perihelion values, T rot = 75 K for perihelion values, and T rot = 60 K for post-perihelion values. Owing to its sensitive dependence on T rot , the mixing ratio for CH 4 is given as a range based on the values in

	Year		1998	2005		2018 d
	Phase	Pre-perihelion	Pre-	Pre-	Perihelion	Post-	Mean among	Range in
				perihelion	perihelion		perihelion	Comets e	Comets f
	CO	< 3.2 a	10±6 b	-	1.72±0.12	-	1.26±0.23	6.1±1.6(19)	0.30 -26
	CH 4 C 2 H 6	-< 0.08 a 0.22±0.13 b 0.14±0.02 c --	0.63 -1.52 0.24±0.04	-0.29±0.02 d	<0.55 -<0.89 0.14±0.03	0.88±0.10(19) 0.55±0.08(27)	0.11 -1.6 0.037 -1.9
	Notes.

a Abundances taken from

Weaver et al. (1999)

.

b Abundances taken from

Mumma et al. (2000)

. c Abundances taken from DiSanti et al. (2013). d This work.

Table 2 .

 2 Weighted

mean values of Q(H 2 O) were (2.86  0.15) x 10 28 mol s -1 pre-perihelion, (3.05  0.

09) x 10 28 mol s -1 near perihelion, and (2.02  0.25) x 10 28 mol s -1 post-perihelion.

Acknowledgements

Data for this study were obtained at the NASA Infrared Telescope Facility (IRTF), operated by

the Infrared Telescope Facility, which is operated by the University of Hawaii under contract NNH14CK55B with the National Aeronautics and Space Administration

), and International Space Science Institute Team 361. We thank an anonymous reviewer for suggestions that improved the paper. We acknowledge and thank the entire staff at IRTF for their support during our observations of G-Z.


Thomas Schmitt
email: thomas.schmitt@ecp.fr

Gabriel Sta↵elbach

Sébastien Ducruix

Stefan Gröning

Justin Hardi

Michael Oschwald

Large-Eddy Simulations of a sub-scale liquid rocket combustor: influence of fuel injection temperature on thermo-acoustic stability

The present work aims at the prediction of high frequency combustion instabilities in cryogenic rocket engines using Large Eddy Simulation (LES). The sub-scale rocket combustor BKD from DLR has 42 coaxial injectors fed with liquid oxygen and gaseous hydrogen operating at conditions representative of full-scale engines. The combustor is found to be stable for low hydrogen injection temperatures and unstable for higher ones. The LES solver is first validated on the Mascotte test-bench from ONERA. Two BKD cases are then simulated, corresponding to naturally stable and unstable regimes in the experiment. The recent progresses and results are described and analyzed.

Introduction

Combustion instability arises from the resonant coupling between the unsteady heat release rate of the flame and the acoustic eigenmodes of the system. The early developments of liquid propellant rocket engines during the twentieth century have often led to large-scale instabilities with serious consequences. The very high power densities in these motors (of the order of 50 GW m 3 ) can induce extreme levels of pressure fluctuations during unstable episodes, which increase heat fluxes towards the wall, damaging the chamber and leading in most cases to the failure of the launch. High-frequency instabilities, corresponding to the coupling with a transverse acoustic mode of the combustion chamber, are the most dangerous ones and remain an issue during the development of liquid rocket engines, because of the curent inability to accurately predict their occurence. The DLR (German Space Center) has developed a cylindrical combustion chamber equipped with 42 coaxial injectors and operating in the transcritical regime (ie at a pressure above the critical pressure of oxygen but an oxygen injection temperature lower than the critical temperature). This sub-scale thrust chamber, with a power up to 91 MW, is referred to as BKD. [START_REF] Gröning | Influence of hydrogen temperature on the acoustics of a rocket engine combustion chamber operated with lox/h2 at representative conditions[END_REF][START_REF] Gröning | Influence of hydrogen temperature on the stability of a rocket engine combustor operated with hydrogen and oxygen[END_REF][START_REF] Gröning | Injector-driven combustion instabilities in a hydrogen/oxygen rocket combustor[END_REF] Since this experiment can operate both under stable and unstable operating conditions, it provides a unique opportunity to assess the ability of the simulation to predict and discriminate between these two situations.

There is, to the authors' knowledge, only one three-dimensional numerical study of a complete engine representative of liquid propellant rocket engines available to the scientific community. These are the calculations made with the AVBP solver. [START_REF] Urbano | Study of flame response to transverse acoustic modes from the les of a 42-injector rocket engine[END_REF][START_REF] Urbano | Exploration of combustion instability triggering using large eddy simulation of a multiple injector liquid rocket engine[END_REF] These first calculations made it possible to show the feasibility of a complete simulation of the BKD bench. The results are particularly remarkable because they allow to identify the most unstable cases (in agreement with the experiments) and the structure of the associated unstable modes. However, the thermo-acoustic instability did not appear naturally in the computation, unlike in the experiments. An initial pulse of variable amplitude was added during the numerical tests to determine the so-called non-linear stability of the system.

The aim of this work is now to predict the natural instability, i.e. without any additional triggering of the system. The reference BKD test case is the one presented in Gröning et al., [START_REF] Gröning | Influence of hydrogen temperature on the acoustics of a rocket engine combustion chamber operated with lox/h2 at representative conditions[END_REF][START_REF] Gröning | Influence of hydrogen temperature on the stability of a rocket engine combustor operated with hydrogen and oxygen[END_REF] where the temperature of the hydrogen flow is gradually increased during the experiment, demonstrating higher pressure oscillations for the higher temperatures (flow rates being maintained). The objective is then to evaluate the ability of high-performance computing to predict the occurrence of combustion instabilities in a rocket engine via a change in operating conditions, an aspect of primary importance for the space industry. The AVBP solver from Cerfacs and IFPEN used here is first briefly presented in Sec. 2. Validations based on Large-Eddy Simulations of the Mascotte test-bench from Onera 26 are then shown in Sec. 3. Finally, two cases from the BKD test bench are simulated and detailed in Sec. 4, corresponding to naturally stable and unstable regimes in the experiment. Conclusions and perspectives are drawn in Sec. 5.

Governing equations, models and numerics

The Favre-filtered, fully compressible Navier-Stokes equations are used to formulate the Large-Eddy Simulations (LES) approach. [START_REF] Poinsot | Theoretical and numerical combustion -Third Edition[END_REF] The species and heat fluxes use classical gradient approaches. The fluid viscosity and the heat di↵usion coe cient are calculated following the Chung et al. method, [START_REF] Chung | Generalized multiparameter correlation for nonpolar and polar fluid transport properties[END_REF] and the species di↵usion coe cients are obtained via a Schmidt number along with the Hirschfelder-Curtiss approximation. [START_REF] Poinsot | Theoretical and numerical combustion -Third Edition[END_REF] The error associated to the constant Schmidt numbers assumption mostly appears in zones of high-density gradient, [START_REF] Oefelein | Mixing and combustion of cryogenic oxygen-hydrogen shear-coaxial jet flames at supercritical pressure[END_REF] where it turns to be negligible in comparison with the dominant turbulent viscosity. In the flame region, no important variation of the Schmidt numbers have been observed. [START_REF] Ribert | Counterflow di↵usion flames of general fluids: Oxygen/hydrogen mixtures[END_REF] Following DNS results of Oefelein, [START_REF] Oefelein | Mixing and combustion of cryogenic oxygen-hydrogen shear-coaxial jet flames at supercritical pressure[END_REF] the Soret and Dufour e↵ects, which are shown to be much smaller than ordinary di↵usion, are neglected.

The subgrid-scale (SGS) stress tensor is modeled with the WALE model, [START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient[END_REF] well-suited for shear flows. The SGS energy and species fluxes are modeled using the gradient transport assumption, introducing turbulent Prandtl and Schmidt numbers, all set to 0.6. In the absence of fully validated models, additional SGS terms, depicted in Taskinoglu et al. [START_REF] Taskinoglu | A posteriori study using a DNS database describing fluid disintegration and binaryspecies mixing under supercritical pressure: heptane and nitrogen[END_REF] and that can be neglected in ideal case, are not considered here. Thermodynamic non-idealities are accounted for using the Soavre-Redlich-Kwong equation of state, [START_REF] Soave | Equilibrium constants from a modified Redlich-Kwong equation of state[END_REF] used for a consistent derivation of all thermodynamic coe cients and functions. [START_REF] Bellan | Theory, modeling and analysis of turbulent supercritical mixing[END_REF][START_REF] Oefelein | Mixing and combustion of cryogenic oxygen-hydrogen shear-coaxial jet flames at supercritical pressure[END_REF] The compressible unstructured solver AVBP [START_REF] Moureau | High-order methods for DNS and LES of compressible multi-component reacting flows on fixed and moving grids[END_REF][START_REF] Schönfeld | Steady and unsteady flows simulations using the hybrid flow solver avbp[END_REF] developed by CERFACS and IFP-EN, with specific modules by EM2C for real gas e↵ects is used for this study. Its Taylor-Galerkin weighted residual central distribution scheme, called TTGC, is third-order in time and space. [START_REF] Colin | Development of high-order Taylor-Galerkin schemes for unsteady calculations[END_REF] The Jacobian matrices of non-viscous fluxes used by the scheme are written to be fully consistent with the non-ideal thermodynamics. Boundaries are treated with the characteristic wave decomposition method NSCBC derived for non-ideal thermodynamics. [START_REF] Okong | Consistent boundary conditions for multicomponent real gas mixtures based on characteristic waves[END_REF] Chemical conversion is handled using the infinitely fast chemistry model. [START_REF] Schmitt | Large-eddy simulation of oxygen/methane flames under transcritical conditions[END_REF] This model is modified to account for the following species: H2, O2 , H2O and OH in order to properly calculate the burnt gas temperature. Equilibrium reference mass fractions are tabulated in terms of the mixture fraction and its variance, both transported during the simulations. The species source terms are then computed following the method described in Schmitt et al. [START_REF] Schmitt | Large-eddy simulation of oxygen/methane flames under transcritical conditions[END_REF] 3. Validation of the LES solver: simulation of the Mascotte test cases

Experimental configurations

The validation configuration is the Mascotte experiment operated by Onera, [START_REF] Vingert | The mascotte single injector 60 bar hot test for code validation -testcase specifications and some more experimental data[END_REF] which has been extensively used for experimental studies of cryogenic combustion. [START_REF] Candel | Structure and dynamics of cryogenic flames at supercritical pressure[END_REF][START_REF] Habiballah | Experimental studies of high-pressure cryogenic flames on the Mascotte facility[END_REF] The present simulations reproduce cases A60 and C60, corresponding to transcritical flames. [START_REF] Candel | Structure and dynamics of cryogenic flames at supercritical pressure[END_REF][START_REF] Juniper | Structure of cryogenic flames at elevated pressures[END_REF] A unique coaxial injector produces a dense oxygen stream (i.e. liquid-like oxygen) at low velocity (less than 10 m/s), surrounded by a high-velocity gaseous hydrogen stream (more than 100 m/s), in a chamber at pressure that is above the critical one for both reactants. For both cases the chamber pressure is set to 60 bar. Cases A60 and C60 di↵er only by the hydrogen mass flow rate, which is 50% larger for case A60. It leads to a larger momentum flux ratio between fuel and oxygen for case A60 compared with case C60. [START_REF] Candel | Structure and dynamics of cryogenic flames at supercritical pressure[END_REF] Oxygen is injected at 83 K, below its critical temperature (T C,O 2 = 154 K), and is transcritical, while hydrogen, injected at 300 K, is supercritical and gaseous (T C,H 2 = 33 K). Under such conditions, the density of oxygen (⇢ O 2 = 1340 kg.m 3 ) is much larger than that of hydrogen (⇢ H 2 = 5.25 kg.m 3 ), leading to a density ratio S ⇢ between the reactants close to 250 at injection. Note that in both cases, hydrogen is always in excess.

Computational domain, mesh and boundary conditions

The computational domain is sketched in Fig. 1a and reproduces the geometry of the experimental chamber. The section of the combustion chamber is square-shaped (50 mm x 50 mm), and the domain length is of the order of 1 000 mm. A zoom close to the coaxial injector is shown in Fig. 1b. The thin (o(0.1 mm)) tapered lips are also accounted for in the simulation. The three-dimensional simulations are performed on a mesh containing 2,000,000 points and 11,000,000 tetrahedra. The smallest cell is located at the tip of the separator between oxygen and hydrogen, meshed with 5 cells in its thickness, which proved to be su cient for a proper flame stabilization (Fig. 1b). The mesh is progressively coarsened downstream, where the flow and flame structures are less sharp. Walls are treated as adiabatic with a wall-law boundary condition. The inlet and outlet conditions are both set with non-reflecting boundary conditions, with a relaxation on the pressure at the exit boundary condition. Anisotropic turbulent velocity fluctuations are superimposed to the bulk flow at the injection (on both oxygen and hydrogen streams), with an intensity of 3 %.

Results and comparison with experimental data

Iso-surfaces of temperature (2000 K) colored by the axial velocity are presented in Fig. 2. It qualitatively shows the flame brush. The flame topology is similar in both cases, with an initially low spreading angle in the near injector region followed by a sudden expansion region (x ⇡ 10d). Finally, large scale structures develop in the second part of the flame as the flow interacts with the lateral walls. The flame in case C60 is much longer than the one obtained in case A60. The lower length in case A60 is expected to come from the much larger momentum flux ratio for case A60 compared with case C60, increasing the shear and the combustion e ciency. [START_REF] Candel | Structure and dynamics of cryogenic flames at supercritical pressure[END_REF] Numerical results are compared with experimental results from Juniper et al. [START_REF] Juniper | Structure of cryogenic flames at elevated pressures[END_REF] in Fig. 3. Iso-contours of OH mass fraction (white lines) from LES are super-imposed with Abel transform of the time-averaged natural OH* emission images (in color). Such a comparison is qualitative; only the flame position can be properly compared. In both cases, the flame shape is properly retrieved and the axial position of maximum spreading position is correct. Note that in both calculations and experiments, the flame is closed at the end of the visualisation windows for case A60 while it is still open for case C60. The overall agreement between simulation and expriment is very good.

Large Eddy Simulation of BKD 4.1 Experimental configurations

The BKD experiment [START_REF] Gröning | Influence of hydrogen temperature on the acoustics of a rocket engine combustion chamber operated with lox/h2 at representative conditions[END_REF][START_REF] Gröning | Influence of hydrogen temperature on the stability of a rocket engine combustor operated with hydrogen and oxygen[END_REF][START_REF] Gröning | Injector-driven combustion instabilities in a hydrogen/oxygen rocket combustor[END_REF] is operated at the European Research and Technology Test Facility P8 for cryogenic rocket engines at DLR Lampoldshausen is a sub-scale thrust chamber featuring 42 H 2 /O 2 coaxial injectors (Fig. 4a). This configuration has the particularity of featuring self-excited combustion instabilities depending on the injection conditions (Fig. 4b). The chamber is instrumented with a high frequency measurement ring including 2 dynamic pressure sensors, located at 5 mm from the injection plan. There is no optical access, but the instrumentation is su cient for the detection of unstable regimes and the modal identification of the observed instabilities. [START_REF] Gröning | Injector-driven combustion instabilities in a hydrogen/oxygen rocket combustor[END_REF] The chamber walls are water-cooled to limit their maximum temperature and heat loss measurements are available. The length and the radius of the chamber are of the order of 100 mm, while coaxial injectors features characteristic diameters representative of In this study, we concentrate on the influence of hydrogen temperature injection on the thermo-acoustic stability of the chamber. [START_REF] Gröning | Influence of hydrogen temperature on the stability of a rocket engine combustor operated with hydrogen and oxygen[END_REF] All the other injection characteristics (hydrogen and oxygen mass flow rates, oxygen temperature) are kept unchanged during the simulations. Injection conditions for the two cases are shown in Tab. 1. The chamber pressure for both cases is close to 80 bar. Oxygen is transcritical (⇢ LOx = 1100 kg m 3 ) and injected at low velocity (O(10 m s 1 )), while hydrogen supercritical (⇢ H 2 = 29 kg m 3 ) and injected at high velocity (O(100 m s 1 )).

The two cases retained for the simulations are represented in Fig. 4b. Case T57 leads to lower pressure fluctuation amplitudes in the chamber than case T131. In particular, for case T57, the second transverse (2T) acoustic mode slightly dominates the first transverse (1T) one. As a consequence, case T57 is considered being stable for the 1T mode and moderately excited for the 2T mode. On the contrary, case T131 presents large amplitudes of fluctuations of the 1T mode (compared with case T57) and the mode is considered as excited, while the 2T mode is stable. In the following, for the sake of simplicity, case T57 is said to be stable and case T131 is considered unstable.

Computational domain, mesh and boundary conditions

The computational domain is shown in Fig. 5 and corresponds to the experimental chamber, including the hydrogen and oxygen domes and feeding lines. In order to properly represent the outlet flow and acoustic impedance, the outlet nozzle is also simulated. The grid (Fig. 6) comprises 55,000,000 points and 300,000,000 tetrahedra. It is highly refined close to the annular injectors that feature very small geometrical elements (0(0.1mm)). All the solid boundaries are treated as adiabatic no-slip walls and iso-thermal (500 K) wall-law are used at the chamber walls to account for water cooling. Oxygen and hydrogen are injected in the domes using non-reflecting boundary conditions with mass flow rates and temperatures imposed. As the nozzle is choked, nothing is imposed at outlet.

Temporal evolutions

The temporal evolution of pressure in the chamber is plotted in Fig. 9 for both cases. Frequencies lower than 5 000 Hz have been removed thanks to a third-order Butterworth filter. These plots show the temporal evolution since the beginning of the simulation on the final grid (after the interpolation from the coarser grid). After a few milliseconds, the pressure in case T57 reaches a quasi-steady behavior, with pressure oscillations around 4 bar peak-to-peak. After 10 ms, the amplitude of the pressure oscillations increases suddenly and the system becomes naturally unstable. After 17 ms, the system seems to have reached a limit cycle, with pressure oscillations around 13 bar peak-to-peak. On the contrary, case T131 remains fully stable during the whole simulation, with pressure oscillations lower than 3 bar peak-to-peak. Figure 10a shows the power spectral density of pressure in the chamber for case T57, once the instability has reached the limit cycle. Two peaks dominate. The first mode at 10 kHz corresponds to the first transverse mode of the chamber. This is shown in the bandpass-filtered signal (around 10 kHz) in Fig. 10b, where the signals from P1 and P3 (see Fig. 11(1)) are in phase opposition and with a much larger ampltiude that the ones from P2 and P4. The second mode at 20 kHz seems to correspond to the first radial mode of the chamber. [START_REF] Urbano | Exploration of combustion instability triggering using large eddy simulation of a multiple injector liquid rocket engine[END_REF] Successive snapshots are shown in Fig. 11. These results are in opposition with those from experiments (T57 should be stable for the 1T mode and T131 unstable). The reason for this discrepancy is under investigation and is expected to come from a insu cient grid resolution in the annular injector. Interestingly, the unstable mode obtained in LES share strong similarities with the unstable one in the experiments 7 : same dominant modes (1T and 1R), similar frequencies, similar pressure fluctuations amplitude. This may be due to the little di↵erences observed between the flow fields for the two cases, as shown in Fig. 8. Still, these results are encouraging, because this is the first time a natural transition from a stable to an unstable behavior is obtained by numerical simulations in conditions representative of LOx/H 2 rocket engines (i.e. a multi-injector configuration under LOx/H 2 transcritical conditions) without external triggering.

Influence of thermo-acoustic instability on flame topology (case T57)

Figure 12 shows instantaneous fields of temperature for case T57 after the growth of the instability (22 ms on Fig. 9a).

Compared with the previous illustrations shown in Sec. 4.3, the flames are now much shorter, in particular in the center of the chamber where the transverse acoustic velocity is maximum. These flames, located in a pressure node and a velocity anti-node, show a sinuous motion. [START_REF] Hakim | Large eddy simulations of multiple transcritical coaxial flames submitted to a high-frequency transverse acoustic modulation[END_REF][START_REF] Hakim | Dynamics of a transcritical coaxial flame under a high-frequency transverse acoustic forcing: Influence of the modulation frequency on the flame response[END_REF] The outer flames, on the contrary, essentially located in a pressure anti-node leading to inlet velocity modulations, [START_REF] Urbano | Study of flame response to transverse acoustic modes from the les of a 42-injector rocket engine[END_REF] present a varicose-like shape.

Conclusions

This paper present work-in-progress calculations of the BKD experiment from DLR. This sub-scale combustor is naturally unstable depending on the injection conditions. In particular, it has been shown 6 that the configuration is stable for low hydrogen temperature injection and become unstable for a higher temperature. Two cases are then retained for the simulations: one stable (T57) and the other one unstable (T131). These cases only di↵er by the hydrogen injection temperature. These simulations are performed in a LES framework with the AVBP RG solver. This solver is first validated simulating the Mascotte configurations from ONERA operating at supercritical pressure. Very good agreements are obtained on the configurations A60 and C60. Current numerical simulations predict that case T57 is unstable and case T131 is stable. Interestingly, the unstable modes obtained in LES share strong similarities with the unstable ones in the experiments [START_REF] Gröning | Injector-driven combustion instabilities in a hydrogen/oxygen rocket combustor[END_REF] : same dominant modes (1T and 1R), similar frequencies, similar pressure fluctuations amplitude. Although not in agreement with experimental measurements, current results are very encouraging as, to the authors knowledge, this is the first time a combustion instability is naturally predicted by a LES solver in a configuration representative of a LOx/H 2 rocket engine (transcritical regime, large number of injectors, cylindrical chamber). Investigation are now focusing on the injection flow, as it may be a key element for the explanation of the observed discrepancy.

Figure 1 :

 1 Figure 1: (a) Sketch of the computational domain. (b) Longitudinal slice of the mesh

Figure 2 :Figure 3 :

 23 Figure 2: Iso-surfaces of temperature (2000 K) colored by axial velocity (blue: -25 m/s ; red: 150 m/s) (a) Case A60, (b) Case C60.

Figure 4 :

 4 Figure 4: (a) Sketch of the BKD chamber 5 (b) Experimental results during ramp tests, 5 where hydrogen temperature is slowly changed. Gliding RMS calculation of the dynamic pressure sensor signals in the chamber (5 mm from injection plan) for the 1T and 2T modes, dynamic pressure sensor signal in chamber (5 mm from injection plan) and injection conditions.

Figure 8 :

 8 Figure 8: Longitudinal slices of instantaneous fields of (a) temperature (blue: 50 K ; red: 3 600 K), (b) axial velocity (blue: -250 m/s ; red: 500 m/s) and (c) density (blue: 5 kg/m 3 ; red: 1150 kg/m 3 ). Top: Case T57. Bottom: Case T131.

Figure 9 :

 9 Figure 9: Pressure evolution in chamber for (a) case T57 and (b) case T131. The signal is band filtered between 5 kHz and 50 kHz using a third order Butterworth filter.

Figure 10 :

 10 Figure 10: Case T57. (a) Power spectral density of the pressure signal in the chamber (Welch's method, with 5 blocks and 50 % overlapping, df=50 Hz) computed after 17 ms. (b) Bandpass-filtered (around 10 kHz) pressure signal (at 18 ms) in the chamber for 4 di↵erent probes positioned along the chamber wall at 5 mm from the injection plan and 90°f rom each other (see Fig. 11(1)).

Figure 11 :

 11 Figure 11: Case T57. Transverse slices (5 mm from injector plan) of instantaneous fields of pressure (Blue: 60 bar ; red: 90 bar) during one period of oscillation of the first mode at 10 kHz.

Figure 12 :

 12 Figure 12: Case T57. Longitudinal slices of instantaneous fields of (a) temperature (blue: 50 K ; red: 3 600 K), (b) axial velocity (blue: -200 m/s ; red: 500 m/s), (c) density (blue: 5 kg/m 3 ; red: 1 150 kg/m 3) and (d) pressure (blue: 60 bar ; red: 90 bar). t = 22 ms.

Table 1 :

 1 Injection conditions for the two selected cases of the BKD experimental setup.[START_REF] Gröning | Influence of hydrogen temperature on the acoustics of a rocket engine combustion chamber operated with lox/h2 at representative conditions[END_REF][START_REF] Gröning | Influence of hydrogen temperature on the stability of a rocket engine combustor operated with hydrogen and oxygen[END_REF][START_REF] Gröning | Injector-driven combustion instabilities in a hydrogen/oxygen rocket combustor[END_REF]

	Case	ṁLOx [kg/s] ṁH 2 [kg/s] T LOx [K] T H 2 [K]	Exp. behaviour
	T57	5.8	0.97	107.6	57.6	Stable 1T mode, moderately excited 2T mode
	T131	5.8	0.97	107.6	131.7	Excited 1T mode, stable 2T mode

Acknowledgements

Support provided by SAFRAN Snecma Space Engines Division, the prime contractor of the Ariane launcher cryogenic propulsion system and CNES, the French National Space Agency, is gratefully acknowledged. This investigation was carried out in the framework of the French-German REST program initiated by CNES, ASL and DLR. Thanks to Dmitry Suslov, Stefan Gröning and the P8 test bench team for performing the experiments. This work was granted access to the HPC resources of IDRIS and CINES made available by GENCI (Grand Equipement National de Calcul Intensif) under the allocation x20162b6176. The support of CINES for access to the computational resources of Occigen under a "grand challenge" allocation is acknowledged. Part of these calculations were performed in 2016 thanks to a special allocation on Curie made available by TGCC.

Both simulation are started on a coarser grid in order to install the large scale recirculation in the chamber. It is then interpolated on the final grid. One convective time is estimated to be close to 5 ms.

Flame tolopogy (stable flames)

Iso-surfaces of temperature (2 100 K) that qualitatively represent the flame brush for case T131 are shown in Fig. 7. The surfaces are highly corrugated by turbulence and extend up to the middle of the chamber. Longitudinal slices of temperature, axial velocity and density are shown for both cases T57 and T131 in Fig. 8. Note that case T57 is here considered before the growth of instability detailed in Sec. 4.4 (at 8 ms on Fig. 9) so that it can be qualitatively compared with case T131. The two cases show a similar flame topology. The cold and dense


T Schmitt
email: thomas.schmitt@ecp.fr

L Hakim

M Boileau

G Staffelbach

A Ruiz

S Ducruix

B Cuenot

S Candel

Large-Eddy Simulation of a multiple injector cryogenic combustor under transcritical conditions and large amplitude high frequency modulations

Keywords: LES, flame dynamics, transcritical combustion, high-frequency combustion instability

.

Introduction

Combustion instabilities have often hindered developments of Liquid Rocket Engines (LRE). Such extreme dynamical phenomena were experienced in many high performance engines. It took a considerable amount of testing to eliminate high frequency instabilities arising in the F-1 of the Saturn V rocket [START_REF] Oefelein | [END_REF]. Such instabilities were also encountered during the Ariane development program leading to the in-flight destruction of the second launcher on November 23, 1980. After a few seconds of operation high-frequency oscillations developed in the first stage propulsion system resulting in a rapid erosion and melting of the injection face due to the rapid intensification of heat fluxes under large amplitude oscillations of the thrust chamber parameters. Because this problem was so pervasive it has been extensively investigated. Most observations were made in real scale engines with relatively few model scale investigations. Theoretical understanding and experimental data based on large scale testing are summarized in [3] while more recent advances are reviewed in [START_REF]Liquid Rocket Engine Combustion Instability[END_REF][START_REF] Culick | Unsteady Motions in Combustion Chambers for Propulsion Systems[END_REF][START_REF] Candel | [END_REF].

The present research project aims at advancing the state of the art in this field by exploring the possibilities of Large-Eddy Simulations (LES) with the following three objectives: (1) Check the numerical methodology by comparing the LES data to the experimental results, (2) Observe the dynamical response of coaxial flames in a high-frequency transverse acoustic field and complement experimental data and (3) Underline differences in flame dynamics depending on whether the flame is submitted to acoustic pressure or acoustic transverse velocity. The simulations reported in what follows correspond to operating conditions prevailing in LREs reproduced in a multiple injector laboratory scale test configuration (MIC) equipped with an external actuator (designated as the VHAM). The VHAM induces large amplitude transverse acoustic modulations which interact with the multiple cryogenic flames established in the MIC. This device is fed with liquid oxygen injected in the central channels of coaxial injectors and this stream is surrounded by gaseous hydrogen or methane. The mean chamber pressure exceeds the critical pressure of the liquid oxygen (p c (O 2 ) = 5.04 MPa) which is injected at a temperature well below its critical value T c (O 2 ) = 154 K. Under these transcritical conditions, the working fluid state departs from that of a perfect gas [7], surface tension vanishes and mass transfer from the dense oxygen stream to the lighter surrounding stream is essentially governed by turbulent production of surface area and mixing. Combustion thus takes place in a complex flow featuring large density gradients and its simulation raises considerable physical and mathematical difficulties. The flame structure under transcritical conditions is now well documented, for example in [8][START_REF] Singla | Proc. Combust. Inst[END_REF][START_REF] Davis | [END_REF][11]. These data provide useful information for numerical modeling and have been used to guide simulation of transcritical flows. In the past decade, efforts have been made to establish computational methodologies for transcritical injection and combustion including real gas thermodynamics, transport models and modified algorithms for numerical integration. LES of turbulent transcritical jets and jet flames have been carried out in the transcritical range by different research groups ( [12][13][14][15][16][START_REF] Ribert | Candel[END_REF][18][START_REF] Matsuyama | 48th AIAA Aerospace Sciences Meeting[END_REF][START_REF] Schmitt | [END_REF]). Much of this work has concerned coaxial flames formed by a transcritical oxygen stream surrounded by a high speed annular flow of hydrogen or methane. These simulations constitute the first step towards the simulations of more complex systems, representative of real engines. The experimental configuration simulated in the present study has served in recent years to obtain an understanding of the driving and coupling processes in high-frequency instabilities. The Multiple-Injector Combustor (MIC), featuring five injectors representative of a LRE injection plane system, was conceived [START_REF] Méry | Combustion for Aerospace Propulsion[END_REF]21], and mounted on the Mascotte cryogenic combustion facility of ONERA. The simulations presented in this article reproduce the MIC equipped with the VHAM actuator that produces transverse acoustic modes in the chamber. The experimental setup is briefly described in Sec. 2. The LES solver and the numerical configuration used in the present investigation are reported in Sec. 3. Results are discussed in Sec. 4. Attention is focused on the flame motion induced by the transverse acoustic modulation.

Brief description of the simulated experimental setup

The experimental setup comprises a rectangular combustor featuring a backplane with five coaxial cryogenic injectors. This system is mounted on the Mascotte test-bench from ONERA [START_REF] Méry | Combustion for Aerospace Propulsion[END_REF].

The system operates at 6.7 MPa and liquid oxygen (LOx) is injected through the inner channel at 80 K. The LOx stream is surrounded by a high velocity flow of gaseous methane (GCH 4 ). The velocities of these LOx and GCH 4 streams are O(1) m.s -1 and O(100) m.s -1 , respectively. Under such conditions, representative of rocket engine steady state conditions of operation, the LOx is transcritical while GCH 4 is supercritical. The operating point conditions are gathered in Tab. 1. The MIC is equipped with an actuator, known as the Very High P ch [MPa] ṁLOx [g.s Amplitude Modulator (VHAM) [START_REF] Méry | Combustion for Aerospace Propulsion[END_REF][START_REF] Méry | Mécanismes d'instabilités de combustion hautefréquence et application aux moteurs-fusées[END_REF], which was designed to generate extreme acoustic levels in the chamber in order to respresent the pressure perturbations arising in combustion instabilities. This system is shown in Fig. 1. It consists of a toothed wheel that alternately blocks the two outlet nozzles of the chamber. The whole mass flow is thus ejected alternately from the top and bottom nozzles. For an appropriate rotation rate of the wheel corresponding to the transverse mode eigenfrequencies, strong resonances are induced in the chamber. Because the nozzles are shut by the wheel in an alternate fashion, transverse modes of odd order are promoted. The MIC is equipped with large quartz windows that give access to the most of the flame region. 3. Numerical setup

Large-Eddy Simulation solver

The unstructured AVBP compressible flow solver (see among others [24,25]) is used in what follows to integrate the threedimensional compressible Navier-Stokes equations for a multicomponent mixture of reactive fluids. The numerical method uses a low-dissipation centered scheme, third-order in space and time [26]. Boundary conditions are treated with the characteristic wave decomposition method NSCBC [27]. A real-gas version of the AVBP flow solver, called AVBP-RG, was recently developed [28,29] with appropriate models to account for real-gas non-idealities. The Soave-Redlich-Kwong equation of state [30] was included in the solver:

p = ρrT 1 -ρb m - ρ 2 a m (T ) 1 + ρb m (1)
with T the temperature and r = R/W, where R designates the universal gas constant and W is the molar mass. The coefficients a m (T ) and b m are calculated according to [START_REF] Poling | The properties of gases and liquids[END_REF]. Equation ( 1) is used for a consistent derivation of the pressure dependence of thermodynamic properties (internal energy, enthalpy, specific heats, compressibility), as described in [START_REF] Poling | The properties of gases and liquids[END_REF][START_REF] Yang | [END_REF]. Chemical conversion is handled using an infinitely fast chemistry model presented in [START_REF] Schmitt | [END_REF]. NSCBC boundary conditions are adapted to real-gas thermodynamics [33]. Because of the steep density and internal energy gradients, real-gas simulations require a specific numerical stabilization procedure described in [29]. The WALE model [34] is used to represent the sub-grid scale (SGS) shear stress tensor. The SGS energy and species fluxes are modeled using gradient transport assumptions with sub-grid Prandtl and Schmidt numbers set to 0.7. Results of calculations in gaseous atmospheric pressure situations as well as in high pressure transcritical cases indicate that the AVBP-RG flow solver in combination with these models provides solutions which are in good agreement with experimental data. These partial validations indicate that one may use this code to examine transcritical combustion dynamics problems.

Numerical configuration

Calculations are carried out in the exact configuration (MIC+VHAM) used in the experimental investigation on the Mascotte test-bench. The mesh shown in Fig. 2 comprises 11,800,000 nodes and 63,300,000 cells. The number of cells used to represent the injector lips where the flame is anchored is equal to 6 which is reasonable in view of previous experience with similar calculations. Walls are treated with an adiabatic no-Figure 2: Full 3D geometry and mesh. ∆ is the local cell size, d no is the nozzle diameter and δ lip and h ch are the lip and chamber widths, respectively. slip condition and an anisotropic turbulence is injected at the inlet following the methodology described in [35]. Bulk and rms profiles, resulting from an independent LES of a long tube with the same Reynolds number and mesh resolution, are used to impose velocity conditions of each injector. The MIC mean pressure is imposed using non-reflecting conditions at the two outlets. The acoustic perturbation generated by the VHAM is represented by an alternate modulation of the two outlets with an imposed pressure:

p UP = A sin(2π f t), (2)
p DOWN = A sin(2π f t + π), (3)
with A the amplitude of the modulation and f its frequency. The present study focuses on the 1T1L modulation (Fig. 3(a)). In this situation, flames are essentially modulated by a transverse acoustic velocity. Acoustic pressure fluctuations are zero for the center flame, and modulate the lateral flames. This approach allows to reach the modulation levels observed in the Mascotte experiments and recorded using dynamic pressure sensors located on the upper and lower side walls (Fig. 3(b)). The resulting acoustic modes suitably match the experimental data, as the frequency of modulation is successively made to coincide with the eigenfrequencies of the system.

Results and discussion

Simulation results are first briefly described and then compared with experimental visualizations in Sec. 4.1. The discussion is then focused on the unsteady flame motion in Sec. 4.2.

Simulation results and comparison to experimental data Flame structure

It is first convenient to examine cuts of the velocity amplitude (in red) together with density (greater than 100 kg.m -3 , in blue) shown in Fig. 4. In the absence of acoustic modulation, the typical transcritical flame geometry is observed [36]. The flame, localized between the oxygen and methane streams, has a limited spreading rate but extends far downstream in the chamber. This flame and the inner jet are strongly sheared by the high velocity gaseous annular flow. The latter quickly transitions to a fully turbulent flow. By opposition, the inner high density jet remains cohesive and penetrates up to 20 diameters into the chamber. Each flame is essentially axisymmetric (Fig. 4(a)) in this case. As acoustic modulation is activated, the jets structure changes in a substantial fashion (Fig. 4(b)). The inner jet length is strongly reduced and alternate shedding of large scale vortices takes place in the annular gaseous stream. The methane jets are violently shaken by the transverse acoustic velocity generated by the modulation. The inner jets (and the flames) no longer remain axisymmetric, but are now flattened in the spanwise direction by the transverse velocity modulation.

Comparison with experimental visualizations

Figure 5 compares LES results to instantaneous visualizations from the experiment. Experimental images correspond to a direct view of the flames, which integrates the light emitted along the line of sight. This is qualitatively reproduced in the numerical results by using a volumetric rendering of the temperature field (in blue). The density field is also highlighted with a volumetric rendering in black. This comparison indicates that simulations are qualitatively able to retrieve what is observed experimentally, in terms of flame length and spreading rate. The two peripheral flames are slightly bent towards the chamber centerline. This is also visible in the experiments (Fig. 5(b)). This feature can be related to the formation of recirculation zones in the upper and lower left corners of the chamber and to a reduced pressure re-gion in the center of the chamber. When transverse acoustic perturbations are generated in the system (Fig. 5(c-d)), flames are shortened in both experiment and simulation. It is also found that the center flames are longer than the outer flames. The lateral flames bending towards the chamber centerplane is noticeably larger when the acoustic modulation is on. In Fig. 5(d), the dense cores become visible (dark region in the center of each flame) and two thick layers of intense luminosity form on their upper and lower sides. This phenomenon can be attributed to the flames flattening (Fig. 4). Both in the LES and in the experiments, the five flames oscillate as a block in the transverse direction following the transverse acoustic velocity. This aspect is investigated in details in Sec. 4.2.

Flames dynamics

Flame dynamics induced by the acoustic forcing is investigated in this section. In what follows, the analysis is based on phase averaged data. The different flames are referred as -2, -1, 0, +1 and +2 according to the corresponding injector transverse position (see. Fig. 6). It is proposed to analyze the dense oxygen jets motion as it is expected to drive the flames dynamics [START_REF] Hakim | Candel[END_REF]. Each individual jet displacement is analyzed by plotting the jet transverse position y i 0 (where i is the flame number) for various phases of the acoustic modulation along the downstream direction. This position is defined as the transverse coordinate of the barycenter of the oxygen mass fraction profile. Results for the center flame (flame 0) are shown in Fig. 6 (left). For x < 10d, one observes nearly sinusoidal oscillations. The inner jet shows little sensitivity to the transverse modulation. This flag-like motion is no longer visible further downstream. Instead, the jet features a up-down 'block' motion. The same behavior is noticed for the lateral flames. This is exemplified with flame +2 in Fig. 6(b). In addition to the transverse motion, lateral flames are also bent towards the center of the domain. This behavior is much more pronounced under acoustic modulation. Interestingly, it seems that this inclination is super-imposed to the transverse motion and does not alter its dynamics.

In an attempt to quantify the transverse flame displacement, spatially averaged position y i 0 (t) -relative to time average position -is computed over one period:

y i 0 (t) = 1 x 2 -x 1 x 2 x 1 y i 0 (x, t) -y i 0 (x) dx (4)
where [x 1 , x 2 ] is the interval of axial averaging and y i 0 (x) is the time average of y i 0 (x, t). A first simple procedure is to perform the average over the whole flame length: x 1 = 0 and x 2 = 25d. The resulting values are plotted over time during one period in Fig. 7 (symbol ×), for flame 0 (a) and flame +2 (b). A harmonic motion, whose shape is similar to the mode shape induced by the acoustic modulation is clearly identifiable. Of interest is also the transverse position of the last part of the jet, the "block motion" region, defined by x 1 = 15d and x 2 = 25d. Results are gathered in Fig. 7 (symbol •) for both flames 0 and +2. As expected from Fig. 6, the amplitude of the transverse motion is larger in the second part of the jet, but differences between the two average intervals are limited. It is now interesting to compare these motions with the analytical model proposed by Mery et al. [START_REF] Méry | Combustion for Aerospace Propulsion[END_REF] to describe the flame movement induced by the transverse modulation. The model named FAME (Flame Acoustic Motion Equation) [START_REF] Méry | Combustion for Aerospace Propulsion[END_REF] has been deduced from experimental observations. It considers each flame as a thin line put set in motion by the transverse acoustic velocity. The flame transverse position in time is then estimated as a function of the transverse acoustic velocity. Under its simplest form, one can assume that the flame motion exactly matches the transverse acoustic motion. We suppose for simplicity a pure transverse mode here. The transverse acoustic velocity fluctuation is then independent of the axial position and is given by:

v (y, t) = V cos        πy i f l L        sin(ωt) ( 5
)
where v is the transverse acoustic velocity generated by the modulation, V is the corresponding amplitude, y i f l is the initial position of flame i (i.e. the coordinate of the center of the injector), L the lateral dimension of the chamber, t the time and ω = 2π f , with f the modulation frequency. Under the assumption that the transverse motion of the flame is small compared to the wavelength of the acoustic mode, the transverse displacement induced

y i 0 (t) = y i f l + V ω cos        πy i f l L        cos(ωt) (6)
One then obtains for flame 0:

y 0 0 (t) = V ω cos(ωt) (7)
and for flame +2:

y +2 0 (t) = y +2 f l + V ω 0.9563 cos(ωt) (8)
Equations 7 and 8 for flames 0 and +2 are plotted in Fig. 7 (continuous line). The model matches remarkably well the numerical results for both flames, indicating that the displacement model proposed by Mery et al. is able to properly account for the transverse motions executed by the flames when they are interacting with an imposed field of transverse acoustic velocity modulations.

Conclusion

This article reports results of Large Eddy Simulations (LES) of a multiple injector combustor mounted on the Mascotte cryogenic combustion test-rig. The MIC features 5 coaxial injectors forming a linear arrangement. This configuration is representative to some extent of real liquid rocket engines, and operates like these engines under transcritical conditions. The MIC is equipped with an acoustic actuator that produces transverse acoustic modulations, representative of a rocket engine instability. It is shown that calculations are able to retrieve experimental observations, both in presence or in absence of acoustic forcing. Under acoustic modulation, the flame lengths are notably reduced and the flames move in the transverse direction at the acoustic modulation frequency. Moreover, flames are flattened in the spanwise direction which nicely matches experimental observations. This phenomena explain the sudden changes in flame topology observed experimentally. The transverse flame displacement is investigated in detail. It is shown that a small amplitude flag-like motion is present in the near injection region. Further downstream the jets move in "block". The flame transverse displacements are compared with the modeled approach suggested in Mery et al. [START_REF] Méry | Combustion for Aerospace Propulsion[END_REF]. In such a configuration the flames are supposed to follow the acoustic modulation field. Good agreement between simulations and the reduced model are obtained. On-going work deals with the analysis of a second simulation at higher frequency that corresponds to the 1T2L acoustic mode of the chamber. The calculations may also be employed to estimate the Rayleigh source terms which arise in the acoustic energy balance allowing a further discussion of the growth or decay of oscillations.

Figure 1 :

 1 Figure 1: Schematic representation of the MIC equipped with the Very High Amplitude Modulator (VHAM)[START_REF] Méry | [END_REF]

Figure 3 :

 3 Figure 3: (a) Experimental (plain) and numerical (dashed) dynamical pressure signals recorded at HF1 and HF2. (b) 1T1L mode structure in the MIC. The five flames are represented by iso-surfaces of temperature (3000 K). The longitudinal slice shows an instantaneous pressure field (blue: 6.6MPa; red: 6.85MPa).

Figure 4 :

 4 Figure 4: LES of the MIC featuring 5 injectors. Instantaneous longitudinal slice of density and velocity magnitude. (a) Without external acoustic modulation. (b) Submitted to a transverse acoustic modulation.

Figure 5 :

 5 Figure 5: Left: Reconstructed emission images deduced from the simulations using a volume rendering of the temperature field (blue) and the density field (black) without modulation (a), with a 1T1L modulation (c) -Right: Experimental instantaneous visualizations of the flame light emission (b) without modulation, (d) with a 1T1L modulation.

Figure 7 :

 7 Figure 7: Theoretical and numerical flame displacement. -From Eq. 6; × Average flame position for x < 25d; • Average flame position for 15d < x < 25d.

Figure 6 :

 6 Figure 6: Longitudinal phase averaged profiles of oxygen jet transverse position for eight phases over one period. Left: flame 0; Right: flame +2. -Under acoustic modulation; ---Without acoustic modulation; • • • Average position under acoustic modulation

Table 1 :

 1 Operating point chosen for the simulation. P ch is the chamber mean pressure, ṁLOx and ṁCH 4 are the mass flow rates of oxygen and methane, respec-

		-1 ]	ṁCH 4 [g.s -1 ]	MR	J
	6.7	106	100	1.1	3.7

tively. MR is the mixture ratio and J the momentum flux ratio. The injection temperatures are kept constant: T LOx = 80K and T CH 4 = 288K.

Acknowledgements

This work is supported by SAFRAN Snecma Space Engines Division, the prime contractor of the Ariane launcher cryogenic propulsion system, CNES and CNRS in the framework of the French-German program on Rocket Engines STability (REST). We acknowledge PRACE for awarding us access to resource Curie-TN based in France at TGCC, under the allocation 2010030367 -RA0367. The authors also thank P. Scouflaire and L. Vingert for their help during the experimental campaign.