Aomar Bourenane
email: aomar.bourenane@safrangroup.com

Matthieu Arzel

Frédéric Guilloud

Alain Thomas

Shuffled Decoding of Serial Concatenated Convolutional Codes

Keywords: Serial Concatenated Convolutional Codes (SCCC), iterative decoding, shuffled decoding, parallel implementations

Shuffled decoding enables to accelerate the extrinsic information exchange during iterative decoding of concatenated codes. It has already been applied to parallel convolutional codes or low-density parity-check codes. In this article, we propose to apply shuffled decoding to serial concatenation convolutional codes. We take advantage of their systematic encoding to propose an efficient shuffled decoding scheme. Compared to a standard iterative decoding scheme, the convergence of our shuffled implementation is obtained within fewer iterations, each one requiring also less time to be completed. This convergence acceleration yields doubling the throughput. We finally show that doubling the throughput comes at a lower cost than doubling the hardware resources, making this shuffled scheme efficient in term of implementation. For instance, the memory usage is 29% more efficient thanks to our proposal than a baseline scheme, which significantly reduces the power consumption of hardware decoders.

I. INTRODUCTION

Serial Concatenated Convolutional Codes (SCCC) [START_REF] Benedetto | Serial concatenation of block and convolutional codes[END_REF] has been proposed by S.Benedetto and G. Montorsi as an alternative scheme to the Parallel Concatenated Convolutional Codes (PCCC) [2] known as turbo-codes. Their error correcting performance is comparable to their parallel counterpart and to LDPC codes [START_REF] Gallager | Low-density parity-check codes[END_REF]- [START_REF] Graell I Amat | Design and performance analysis of a new class of rate compatible serially concatenated convolutional codes[END_REF], while exhibiting a better error floor performance. Moreover, the low complexity of SCCC encoders make them suitable for on-board satellite implementations [START_REF] Benedetto | Mhoms: high-speed acm modem for satellite applications[END_REF] [START_REF]BOOK Recommended Standard CCSDS[END_REF]. SCCC encoding consists in the serial concatenation of two convolutional encoders linked by an interleaver, the first one corresponding to the outer code and the second one to the inner code [START_REF] Benedetto | Serial concatenation of block and convolutional codes[END_REF]. Thus the iterative decoding scheme consists in decoding successively the inner and the outer codewords using Soft-Input Soft-Output (SISO) algorithms and exchanging extrinsic information with each other.

In the literature, several works have been carried out in order to achieve high throughput turbo-code decoders using different parallelism techniques. Most of the contributions concern PC-CCs [START_REF] Weithoffer | 25 Years of Turbo Codes: From Mb/s to beyond 100 Gb/s[END_REF] and few are dedicated to SCCCs decoders even though their efficient implementation is an issue for satellite telemetry systems [START_REF] Bertolucci | A comprehensive trade-off analysis on the ccsds 131.2-b-1 extended modcod (sccc-x) implementation[END_REF]. In [START_REF] Martina | High throughput implementation of an adaptive serial concatenation turbo decoder[END_REF], a parallel SCCC-decoder architecture with 16 concurrent SISO decoders was proposed based on sliding windows parallelism technique [START_REF] Schurgers | Memory optimization of map turbo decoder algorithms[END_REF]. In [START_REF] Shoup | Hardware implementation of a high-throughput 64-ppm serial concatenated turbo decoder[END_REF], the hardware implementation of a high throughput Max-Log-Map SCCC turbo-decoder for an optical channel was proposed. Parallelism was introduced both in the state metrics computation within the SISO decoders and by using sliding windows. In these previous contributions, the parallelism introduced in SCCC decoders was intended to speed up the metric calculations but not the extrinsic information exchange which is crucial to further increase the decoder throughput.

In [START_REF] Argon | A parallel decoder for low latency decoding of turbo product codes[END_REF], the authors present a parallel Turbo Product Code (TPC) decoding approach to increase the throughput while keeping the same error performance level. This approach has been generalized on LDPCs and PCCCs under the name of Shuffled decoding technique in [START_REF] Zhang | Shuffled iterative decoding[END_REF]. Shuffled decoding consists in accelerating the extrinsic information exchanges between concurrent decoders operating on the same codeword. More precisely, the aim is to exchange the extrinsic information as soon as it is produced and not once the component code decoding is over. The idea is to speed up the convergence of the iterative decoding, and hence the throughput of the decoder. However, care has to be taken to design the decoder since increasing the pace of extrinsic information exchange might also decrease the error rate performance. In [START_REF] Lu | A parallel decoder design for low latency turbo decoding[END_REF], the authors presents a parallel decoding design for a PCCC decoder based on the work presented in [START_REF] Argon | A parallel decoder for low latency decoding of turbo product codes[END_REF]. In [START_REF] Muller | Exploring parallel processing levels for convolutional turbo decoding[END_REF], the authors explore and analyze parallelism techniques in parallel turbo decoding, including the shuffled decoding technique.

To the best of our knowledge, shuffled decoding for SCCC has never been addressed in the literature. In this paper, we propose a novel SCCC-based shuffled decoding scheme that increases the processing throughput while improving hardware efficiency.

The rest of the paper is organized as follows. In Section II, we introduce the notations used throughout the paper and we describe the encoding and decoding schemes of SCCC. In Section III, we propose an efficient shuffled implementation for systematic SCCC. In Section IV, the performance of the proposed scheme is addressed and compared to that of a baseline implementation. Conclusions and perspectives are discussed in Section V.

II. MODEL AND NOTATIONS

For the sake of simplicity, throughout the paper we shall refer to the communication channel shown in Fig. 1 ) is alternatively equal to the systematic bit (resp. ) and the redundant bit (resp. ). The SCCC decoder is depicted in Fig. 2. It consists of two serially concatenated SISO decoders denoted outer SISO and inner SISO, an interleaver and a deinterleaver. The trellis length of the outer (resp. inner) code is (resp. 2 ) sections. LLRs can be considered as information regarding the decision to make when estimating a bit. Hereafter, information regarding a bit will be denoted ( ) and decomposed according to: ( ) = ( ) + ( ) + ( ), where denotes the channel information, denotes the a priori information and denotes the extrinsic information. The inner SISO is fed by the channel LLRs ( ). Since the inner encoder is systematic, they correspond alternatively to the inner code systematic bits ( ) and to the inner code redundant bits ( ). The inner SISO is also fed by a priori information about the inner code systematic bits, coming from the outer SISO. Only information related to the systematic bits ( ) is output to feed through the de-interleaver the outer SISO since they correspond to the coded bits of the outer code. The outer SISO computes information on both the systematic bits ( ) and the redundant bits ( ) to feed back extrinsic information to the inner SISO. ( ) is used at the end of the iterative decoding to estimate the transmitted bits. As illustrated in Fig. 2 by the blue dashed wires, we emphasize that part of the channel information, ( ), is input to both the inner and the outer decoders, since the inner code is systematic.

III. SHUFFLED SCCC DECODER IMPLEMENTATION

A. Baseline decoder

A standard implementation for SCCC decoding is illustrated in Fig. 3. In this implementation, the hardware resources are shared between the inner and the outer decoders: a single SISO component is instantiated, as well as a single memory for the state metrics and a single memory for the extrinsic information. These resources are used alternatively to process the inner and the outer decoding, one trellis section being processed every clock cycle. The control unit schedules the computing and drives the memory access depending on the component code being decoded. These resources have thus to be sized to support the greedier case.

As computational resources are concerned, the outer SISO has to calculate 2 pieces of extrinsic information per trellis section: one related to systematic bits ( ) and one related to redundant bits ( ). On the other hand, the inner SISO has to calculate only one extrinsic information which is related to systematic bits ( ): it means that only half of the resources are then required to calculate the extrinsic information related to the inner decoder.

As memory resources are concerned, since the trellis length of the inner code is twice the one of the outer code (the code rate is set to 1/2), so has to be the ( , ) State Metrics memory of the inner decoder compared to the one of the outer decoder. Therefore, in this baseline decoder, half of this memory is not used during the outer decoder processing.

B. Asymmetrical Shuffled Decoding (ASD)

The idea of a shuffled implementation is to speed up the extrinsic information exchange, if possible as soon as it is calculated. In the baseline decoder, extrinsic information is used in one component code decoder after it has been calculated for all the bits involved in the other code component decoder. In the proposed shuffled implementation, both component code decoders will proceed simultaneously, using the extrinsic information produced in parallel by each decoder as soon as it is output. Hence, two SISO decoders have to be instantiated as illustrated in Fig. 4.

The serial concatenation induces that the component decoders do not have the same trellis length to decode (hence the Asymmetrical term in ASD): considering an outer code rate = 1/2, the trellis length of the inner code is twice the outer one's. So processing the component codes does not last the same duration. Let denote the amount of time required to perform a complete inner code decoding. Then, the outer code is decoded within an amount of time equal to /2 since the trellis length is halved.

As illustrated in Fig. 5, we propose to take advantage of the trellis length ratio to run the outer decoder twice while running a single inner decoder. Note that a replica butterfly scheme [START_REF] Zhang | Replica shuffled iterative decoding[END_REF] is implemented to increase the extrinsic information exchange speed and the decoding throughput (as in [START_REF] Weithoffer | 25 Years of Turbo Codes: From Mb/s to beyond 100 Gb/s[END_REF], [START_REF] Weithoffer | On the applicability of trellis compression to turbo-code decoder hardware architectures[END_REF]). In fact, during the second execution of the outer decoder, it generates new extrinsic information ( ) that allows to update up to 50% of the inner decoder a priori information ( ) which can be used in the same iteration. Also, new inner extrinsic information ( ) will provide a more reliable a priori information ( ) to be used by the outer decoder in its second execution. Hence, a complete iteration of the proposed ASD requires an amount of time , when the baseline implementation with a SISO butterfly scheduling runs a complete iteration within an amount of time equal to 3 /2: for the inner code and /2 for the outer one.

Inner SISO ( ) Π ( ) ( ) ( ) ( ) ( ) ˆ Π -1 Π -1 ( ) ( ) ( ) ( ) ( ) ( ) ( ) Outer SISO ( ) DEMUX MUX MUX ( ) ( )

C. Hardware Resources

The decoding of the inner and outer codes will be performed using the Max-log-MAP algorithm [START_REF] Robertson | A comparison of optimal and sub-optimal map decoding algorithms operating in the log domain[END_REF] so as to reduce the computational complexity. To simplify the resource estimation both in terms of memory and of computing units, we assume hereafter that the LLRs and the forward-backward metrics ( , ) are quantized with the same bit-width. The hardware and time resources required by the baseline and proposed ASD schemes are compared in Table I and are detailed hereafter.

1) Memory Requirements: The amount of memory required to save each state metric is given by the number of states per trellis section, i.e. 2 -1 , times the number of sections plus one, that is 2 + 1 for the inner code and + 1 for the outer code. So the number of memory words to store the state metrics of an ASD is given by 2 × 2 -1 (2 + 1) + 2 -1 ( + 1) = 2 (3 + 2). For the baseline scheme, the largest trellis is the inner one and requires 2 (2 +1) memory words to be saved.

The amount of memory required to save the extrinsic information in the inner decoder is given by the number of systematic bits which is 2 . In the outer decoder, it is given by total number of bits since the extrinsic information of both the systematic and the redundant bits have to be saved, that is 2 . Note that the extrinsic information memory requirements are duplicated for the ASD scheme to enable the memory to be simultaneously accessed by the inner and the outer decoders as shown in Fig. 4. In the baseline case, only one extrinsic memory of 2 words is required.

2) Computing Resources: We consider here the amount of resources required to perform the calculation of metrics ( , ) and the extrinsic information in a single trellis section.

The calculation of the state metrics in a trellis section of 2 -1 states is made possible by the use of 2 -1 Add-Compare-Select (ACS) units for or state metrics, that is 2 ACS, both for the inner or the outer decoders since we assume the same constraint length for the two component codes.

The extrinsic information requires first to add the two state metrics to the branch metrics, which requires 2 ADD (Addition operator) per branch, and so 2 × 2 = 2 +1 ADD per trellis section. Then, in the case of the inner decoder (extrinsic information on the systematic bits), the maximum sum over the branches associated to a systematic bit equal to 0 and also over the branches associated to a systematic bit equal to 1 has to be selected. Finding the maximum over 2 /2 values requires 2 -1 -1 Compare-and-Select (CS) units, hence a total of 2 -2 CS. Finally the subtraction between the two previous maximums and subtracting the a priori information and the channel information from the result (cost of 3 ADD) will provide the extrinsic information for the systematic bit on the considered trellis section of the inner decoder. To sum up, the inner trellis calculation requires 2 ACS, 2 +1 + 3 ADD and 2 -2 CS. Grouping the CS and ADD operator into a single ACS operator results finally in a total of 2 +1 -2 ACS and 2 + 5 ADD. As for the outer decoder, the calculation is the same, except that 2 extrinsic information values have to be calculated for each trellis section: one for the systematic bit and one for the redundant bit. Hence the number of operations dedicated to the inner extrinsic information is doubled, yielding a total of 2 + 2 × (2 -2)ACS and 2 × (2 + 5)ADD. Finally the overall complexity of the proposed ASD scheme is obtained by adding the complexity of the inner and outer decoders, resulting in 5 × 2 -6 ACS and 3 × 2 + 15 ADD, whereas for the baseline scheme, the complexity is the one of the outer decoder (3 × 2 -4 ACS and 2 +1 + 10 ADD).

IV. PERFORMANCE COMPARISON

In Fig. 6 both schemes are compared in term of the Bit Error Rate (BER) as a function of signal-to-noise ratio for several number of iterations. The component code is the one described in [START_REF] Benedetto | Mhoms: high-speed acm modem for satellite applications[END_REF] where = 3 (4-state trellis) and = 4320. It is also illustrated in Fig. 7 as a function of the number of iterations for several signal-to-noise ratios. It can be noticed that the proposed ASD scheme converges in less iterations than the baseline scheme. For instance, a BER of 10 -5 is reached at / 0 1.41dB in 8 iterations with the baseline scheme and in only 6 iterations with the proposed ASD scheme. Taking also into account that one iteration of the proposed ASD scheme lasts an amount of time whereas one iteration of the baseline lasts 3 /2, the total decoding time of each scheme is given by ASD = 6 × for the proposed ASD scheme and BD = 8 × 3 /2 = 12 for the baseline decoding scheme. So, the ASD scheme offers a throughput twice that of the baseline decoding scheme. This faster convergence and higher throughput is obtained at the cost of more hardware resources than implemented for the baseline scheme as summed up in Table II with = 3 and = 4320. So the question is to state if this increase in resources is worth the gain in throughput. To this aim, we suggest to calculate the processing efficiency of the proposed architecture when compared to the baseline. For each type of resource, is defined as the ratio between the processing cost (amount of resources times the amount of time they are used) of decoding one codeword with the baseline decoding scheme and that with the ASD scheme:

= BD × BD ASD × ASD (1)
where ASD (resp. BD ) is the amount of hardware resources of the ASD scheme (resp. of the baseline decoding scheme), that is the number of memory words, the number of ADD or the number of ACS used. According to Equation (1) and Table II, the efficiency of the ASD scheme is • 129% in term of memory resources,

• 118% in term of ACS units,

• 133% in term of remaining ADD units. In other words, when compared to the baseline scheme for a targeted throughput, the proposed ASD scheme offers considerable savings in terms of memory and computing resources.

V. CONCLUSION

In this paper, we have proposed a shuffled scheduling of serial concatenated convolutional codes called Asymmetrical Shuffled Decoding (ASD). Shuffled scheduling makes exchange of extrinsic information faster and our results show a significant reduction in the number of iterations required to achieve a given bit error rate. Moreover, the shuffled scheduling enables a parallel implementation of the inner and outer convolutional code decoders. As a consequence, one iteration of our proposed ASD scheme is faster than one of a conventional iterative decoding scheme. We show that the overall throughput is doubled in our ASD implementation at a lower cost than doubling the hardware resources, making ASD efficient in term of implementation. For instance, the memory usage is 29% more efficient thanks to our proposal than the baseline scheme, which significantly reduces the energy consumption of hardware decoders. The generalization of ASD scheme to any code rates will be considered in a future work.

Fig. 1 .

 1 Fig. 1. SCCC encoder and transmission over a binary input AWGN (Bi-AWGN) channel model.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Considered SCCC iterative decoder, where information (LLR) related to the channel observation is coloured in blue.

Fig. 4 .

 4 Fig. 4. Proposed ASD scheme.

Fig. 5 .

 5 Fig. 5. BCJR computation with replica butterfly scheme : (a) inner decoder scheduling, (b) outer decoder scheduling.

Fig. 6 .

 6 Fig. 6. BER performance comparison as a function of / 0 .

Fig. 7 .

 7 Fig. 7. BER performance comparison as a function of number of

TABLE I HARDWARE

 I AND TIME CONSUMPTION OF THE BASELINE DECODING SCHEME VS THE ASD SCHEME

	Criteria		Baseline scheme	Proposed ASD scheme
	SISO processing time	3 /2
	# Memory words	2 (2 +1) +2	2 (3 + 2) + 4
	# ACS		3 × 2 -4	5 × 2 -6
	# ADD		2 +1 + 10	3 × 2 + 15
	10	-5
		0.5	0.65	0.8	0.95	1.1	1.25	1.4
				Eb/N 0 (dB)


J Reynaud
email: jolan.reynaud@onera.fr

P-E Weiss

S Deck

Numerical workflow for scale-resolving computations of space launcher afterbody flows with and without jets

Keywords: ZDES, Hybrid scheme, base flow, jet interactions, compressible flows, AUSM, turbulent flows, space launcher

Scale resolving numerical methods are necessary to accurately predict the afterbody flows of space launchers.

Numerical workflows built for such computations have to be designed in order to obtain a minimal numerical dissipation to resolve fine turbulent structures, a sufficient robustness to capture eventual shock waves and both a reasonable computational cost and an acceptable user workload enabling comparative design studies. This article presents the development of a hybrid numerical framework based on Ducros's sensor and designed to switch from a low-dissipation formulation in presence of vortical structures to a robust formulation around high gradients.

This hybrid workflow is used with ZDES, including its latest automatic mode (ZDES mode 2 (2020)), to simulate a transonic space launcher afterbody experiment with and without a cold air propulsive jet [1] [2]. The salient physical properties of the base flow are investigated and the evaluation of computed results follows the extended nomenclature for validation of simulation techniques [3] from level 0 to 5. The framework is thus validated for instantaneous and mean flow visualizations, base and extension mean pressure coefficient distribution, pressure fluctuation levels, one-point and two-point spectral analyses. Such an efficient automatic RANS/LES strategy could be well suited to study realistic launcher afterbody geometries with propulsive jet(s).

Introduction

Space launcher base flow phenomenology

The prediction of the afterbody flow of a space launcher during its flight is of primary importance for the correct assessment of the aerodynamic performance and the stability of the vehicle, especially in order to establish new designs such as reusable launchers. In flight, the massively separated turbulent flow behind the base and around the nozzle(s) induces two major phenomena: base drag and side-loads. First, base drag, mainly due to the low pressure recirculation zone forming behind the base constitutes a significant contribution to the overall drag. Then, side-loads (related to the buffeting phenomena), caused by pressure fluctuations linked to the large-scale unsteady motion of the recirculation bubble, influence the flight stability and may cause vibrations of the mechanical structures. These adverse effects, especially strong in the high subsonic regime (M ∞ =0.5-0.9), may furthermore be enhanced by the interaction between the propulsive jet(s) and the base flow. This interaction is composed of a blocking effect (nonviscous) and of an entrainment effect (viscous and turbulent). The blocking effect is related to the shape of the jet which acts as a solid wall deflecting the external flow. The entrainment effect tends to empty the base region and thus to decrease base pressure as the jet entrains the external flow through a mixing layer. The main features of the flow are shown on Fig. 1 for a subsonic case of a cylindrical extension behind an axisymmetric forebody.

Review of quantitative studies

The phenomenology qualitatively described in sub-section 1.1 has been studied for decades [START_REF] Delery | Base flows behind missiles[END_REF] and gives a practical framework to study base flows. Howewer it has to be noticed that a thorough quantitative characterization of a space launcher afterbody aerodynamics, for instance according to the extended nomenclature for validation of simulation techniques [START_REF] Weiss | On the coupling of a zonal body-fitted/immersed boundary method with ZDES: Application to the interactions on a realistic space launcher afterbody flow[END_REF] shown in Tab 1, demands high fidelity spatio-temporal methods which are complex and costly to use for both experimental and numerical studies. Hence, the literature showing qualitative results for unsteady quantities of the high subsonic flow around an axisymmetric backward facing step is quite recent as shown in Tab.2 which is not exhaustive but provides an overview of the state of the art. Examining this review, one can see that most efforts have been focused on long extension ratio L D ≥ 1.2 (D: base diameter, L: extension length) representative of Ariane 5 whereas designs presented for the Ariane 6 and Ariane Next [START_REF] De Mirand | Ariane Next, a vision for the next generation of Ariane Launchers[END_REF] rockets indicate choices in favour of much shorter extensions L D ≤ 0.6. These studies also tend to focus on side loads and the base drag issue was given less attention. Furthermore, numerical studies including a propulsive jet are still rare [START_REF] Deck | Unsteadiness of an axisymmetric separating-reattaching flow: Numerical investigation[END_REF] [START_REF] Statnikov | Reduced-order analysis of buffet flow of space launchers[END_REF]. We propose here a numerical study on a relatively short extension geometry (L/D=0.6) with and without a propulsive jet and with interest on both static and dynamic base pressure properties.

Grade Levels of validation 0

Instantaneous flowfield visualization (coherent structure criteria, numerical schlieren)

Integral forces (lift, drag and pitch) 2 Mean aerodynamic field (velocity or pressure profiles)

Second-order statistics (r.m.s. quantities)

One-point spectral analysis (power spectral densities) 5 Two-point spectral analysis (correlation, coherence and phase spectra)

High-order and time-frequency analysis (time-frequency, bicoherence spectra)

Table 1: Extended nomenclature for levels of validation of simulation techniques [START_REF] Weiss | On the coupling of a zonal body-fitted/immersed boundary method with ZDES: Application to the interactions on a realistic space launcher afterbody flow[END_REF] adapted from [START_REF] Sagaut | Large Eddy Simulation for Aerodynamics: Status and perspectives[END_REF].

Study

Physical characterization of high-subsonic base flow

Several physical observations characterizing high-subsonic base flows based on the aforementioned literature have been gathered in this section.

Beginning with static observations, the mean base pressure coefficient (Cp b ) was shown to be mostly radially invariant and does not depend on M ∞ for axisymmetric cases without jet. Values of Cp b are usually measured around -0.15 for case with L/D=1.2 [START_REF] Deprés | Analysis of unsteadiness in afterbody transonic flows[END_REF] [24] an -0.12 for a canonic axisymmetric base flow (L=0) [START_REF] Delery | Base flows behind missiles[END_REF]. Such values testify of an important base drag. In the presence of a propulsive jet, Cp b is lowered [START_REF] Deprés | Analysis of unsteadiness in afterbody transonic flows[END_REF] [10], this decrease is less important when the extension is longer [START_REF] Lê | Etude expérimentale du couplage entre l'écoulement transsonique d'arrière-corps et les charges latérales dans les tuyères propulsives[END_REF]. Regarding the reattachment length of the recirculation bubble, a growth is observed together with M ∞ for cases without jet [22] [24]. For M ∞ =0.8, the reattachment of the recirculation bubble is solid (i.e. occurs on the extension) for L/D≥0.9 and is hybrid or fluid (i.e. occurs downstream of the extension) for shorter extensions [9] [24]. Furthermore, the reattachment length is shortened in the presence of a propulsive jet [22] [24].

Based on spectral analysis and modal decomposition applied on wall pressure signals or numerical data, the dynamic behaviour of an axisymmetric base flow can be decomposed into different main motions. The Strouhal number based on the base diameter St D = f D U∞ is used to characterize fluctuations occurring at frequency f. Authors have observed a low characteristic frequency (St D ≤ 0.1) corresponding to a lateral pumping motion and linked to the growth and decay of the whole recirculation bubble. This frequency is dominant for pressure fluctuations near the base [12][7]. An antisymmetric vertical motion due to a vortex shedding phenomena at a reduced frequency St D ≈ 0.2 has also been found [START_REF] Deprés | Analysis of unsteadiness in afterbody transonic flows[END_REF][6] [START_REF] Weiss | On the dynamics of axisymmetric turbulent separating/reattaching flows[END_REF][7] [START_REF] Van Gent | From particles to pressure PIV-based pressure reconstruction for base flows[END_REF]. The vortex shedding frequency is dominant along the extension and the energy of global effort fluctuations is concentrated near it [START_REF] Weiss | On the dynamics of axisymmetric turbulent separating/reattaching flows[END_REF] [START_REF] Pain | Large scale dynamics of a high Reynolds number axisymmetric separating/reattaching flow[END_REF]. The harmonics of these two low frequency peaks can also be observed on pressure signals. The recirculation region is also moved by fluctuations with a broadband spectrum centered in St D ≈ 0.5-0.6. These fluctuations are attributed to vortical structures formed by Kelvin-Helmoltz instabilities in the shear layer of the recirculation zone [START_REF] Deck | Unsteadiness of an axisymmetric separating-reattaching flow: Numerical investigation[END_REF] [START_REF] Weiss | On the dynamics of axisymmetric turbulent separating/reattaching flows[END_REF]. To complete this description, authors showed that, in the presence of a long extension or a jet, the base flow tends to be stabilized and the amplitude of the antisymmetric oscillations at the vortex shedding frequency is reduced [START_REF] Deprés | Analysis of unsteadiness in afterbody transonic flows[END_REF] [6] [START_REF] Saile | Experimental analysis on near-wake flows of space transportation systems[END_REF]. Moreover, it has been observed that, even though the jet modifies significantly the turbulent kinetic energy distribution of the base flow only for short extensions [START_REF] Van Gent | From particles to pressure PIV-based pressure reconstruction for base flows[END_REF], pressure fluctuations (estimated with Cp rms ) on the base and on the extension are usually higher in the presence of a jet due to its aeroacoustic noise [START_REF] Saile | Experimental analysis on near-wake flows of space transportation systems[END_REF] for all extension lengths.

Turbulence modelling strategies

Thanks to the numerical studies presented in Tab.2 and others on blunt base flows [START_REF] Simon | Zonal-Detached-Eddy Simulation of Projectiles in the subsonic and transonic regimes[END_REF] [START_REF] Simon | Reynolds-Averaged Navier-Stokes/Large-Eddy Simulations of supersonic base flow[END_REF], several results regarding the choice of turbulence modelling strategies for base flow computations can be summarized here.

First, most RANS turbulence models are not well suited for base flow computations. They usually produce an acceptable qualitative flow topology but fail to provide correct mean wall pressure values, recirculation lengths and velocity distribution in the recirculation zone [START_REF] Simon | Zonal-Detached-Eddy Simulation of Projectiles in the subsonic and transonic regimes[END_REF][10] [START_REF] Reynaud | Towards binary gas mixture ZDES for space launcher base flow prediction[END_REF]. The resolution of large eddies driving the dynamics of separated flows is required to increase the representativeness of unsteady base flow simulations [START_REF] Deck | Detached and large eddy simulation of unsteady side-loads over an axisymmetric afterbody[END_REF]. Nevertheless, full Large Eddy Simulations (LES) have almost not been used on such geometries because of the computational cost linked to very fine mesh requirements for attached boundary layers and the long-time sampling needed to capture low frequency dynamics [START_REF] Deck | Detached and large eddy simulation of unsteady side-loads over an axisymmetric afterbody[END_REF]. Strategies enabling to treat attached boundary layer in RANS mode and the separated flow in LES-like mode have then been favoured in the literature [START_REF] Deck | Unsteadiness of an axisymmetric separating-reattaching flow: Numerical investigation[END_REF][21] [START_REF] Pont | Assessment of automatic hybrid RANS/LES models for industrial CFD[END_REF] [23] [START_REF] Scharnowski | Combined experimental and numerical investigation of a transonic space launcher wake[END_REF]. Two main approaches are available for this purpose, providing a good compromise between accuracy and computationnal cost. A first category encompasses hybrid methods in which the framework ensure a continuous transition from RANS to LES equations. Good results on canonical and realistic geometries have been obtained with the ZDES (Zonal Detached Eddy Simulation) approach [START_REF] Deck | Recent improvements in the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF] which is based on the Spalart-Allmaras (SA) model [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF] and has been used for space launcher aerodynamics in [START_REF] Deck | Unsteadiness of an axisymmetric separating-reattaching flow: Numerical investigation[END_REF] [START_REF] Weiss | On the dynamics of axisymmetric turbulent separating/reattaching flows[END_REF][17][3] [START_REF] Hannemann | Launch vehicle base buffeting-recent experimental and numerical investigations[END_REF], with the IDDES (Improved Delayed Detached Eddy Simulation) and the SAS (Scale Adaptive Simulation) approach which are based on Menter's k -ω model [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF] and have been respectively used in [START_REF] Horchler | Hybrid RANS-LES study of transonic flow in the wake of a generic space launch vehicle[END_REF] and [START_REF] Pont | Assessment of automatic hybrid RANS/LES models for industrial CFD[END_REF] and with the PG (Perot and Gadebusch [START_REF] Perot | A self-adapting turbulence model for flow simulation at any mesh resolution[END_REF]) turbulence model which is based on the standard k -and has been used in [START_REF] Pont | Assessment of automatic hybrid RANS/LES models for industrial CFD[END_REF]. Then, a second family of approaches gathers coupling methods where information taken from a RANS zone is used as inlet conditions for a LES zone, often through an overlapping zone.

Such approaches have been successfully used on semi-complex geometries [START_REF] Statnikov | Analysis of pressure perturbation sources on a generic space launcher after-body in supersonic flow using zonal turbulence modeling and dynamic mode decomposition[END_REF][20] [START_REF] Statnikov | Reduced-order analysis of buffet flow of space launchers[END_REF] with a MILES [START_REF] Boris | New insights into large eddy simulation[END_REF] formulation .

The present study will use the ZDES approach for turbulence modelling, including its latest automatic mode (ZDES mode 2 (2020) [START_REF] Deck | Towards an enhanced protection of attached boundary layers in hybrid RANS/LES methods[END_REF]) in order to promote a low user workload numerical framework.

Numerical strategies

For the RANS/LES formulations previously described, the accuracy of the numerical scheme has been shown to be of primary importance. Indeed, reducing the numerical dissipation enables to resolve smaller turbulent scales, to limit the influence of the modelled turbulence diffusion on the computed flow and to avoid the damping of instabilities thus improving the transition from RANS to fully developed LES zones. However, lowering the numerical dissipation tends to decrease the robustness of a method and increases the risk of numerical divergence, especially when strong shocks are present. For this reason, computations including a propulsive jet (containing shock cells) have often used a robust scheme, degrading the resolution of turbulent structures compared to plumeoff cases. Among the numerical studies of Tab.2, for cases without propulsive jet, modified versions of the AUSM+ scheme [START_REF] Liou | A sequel to AUSM: AUSM+[END_REF] including a centered contribution have been used on structured grids with second-order [START_REF] Statnikov | Analysis of pressure perturbation sources on a generic space launcher after-body in supersonic flow using zonal turbulence modeling and dynamic mode decomposition[END_REF] and thirdorder (MUSCL reconstruction) accuracies [START_REF] Weiss | On the dynamics of axisymmetric turbulent separating/reattaching flows[END_REF] [START_REF] Pain | Zonal Detached Eddy Simulation of theFlow Around a Simplified Launcher Afterbody[END_REF]. For cases with jet, for which shock cells have to be resolved, authors have either kept this type of scheme and used very fine grids [START_REF] Statnikov | Reduced-order analysis of buffet flow of space launchers[END_REF](500M points) or switched to more robust schemes like Roe's scheme [10][6] for coarser grids (≈10M points) with a noticeable loss in scale resolution in comparison with jet-off cases. On unstructured grids, Godunov-type methods, even enhanced with high order kexact reconstruction, appeared too dissipative and required very fine grids to solve turbulent structures accurately.

That is why an hybrid centered/upwind version of these methods is proposed in [START_REF] Pont | Assessment of automatic hybrid RANS/LES models for industrial CFD[END_REF]. This scheme uses Ducros's sensor [START_REF] Ducros | Large-eddy simulation of the shock/turbulence interaction[END_REF] to detect vortical regions to apply the centered version of the scheme. A sensor is thus used here to locally reduce numerical dissipation.

As described in [START_REF] Zhao | A general framework for the evaluation of shock-capturing schemes[END_REF], this idea of optimizing dynamically the numerical resolution of the flow based on local properties, especially for turbulent flows containing shocks, is still an open subject of the literature. To develop such a framework, adapted to complex flows on curvilinear geometries, one can study fundamental works on academic shock/turbulence configurations dedicated to the development and the evaluation of high-order shock-capturing formulations [START_REF] Ducros | Large-eddy simulation of the shock/turbulence interaction[END_REF] [START_REF] Garnier | On the use of shock-capturing schemes for large-eddy simulation[END_REF][40] [START_REF] Bogey | A shock-capturing methodology based on adaptative spatial filtering for high-order non-linear computations[END_REF][42] [START_REF] Zhao | A general framework for the evaluation of shock-capturing schemes[END_REF]. The proposed methods are mainly suited for uniform cartesian grids but they provide a thorough analysis of numerical challenges to overcome and promote inspiring sensor formulations that can be customized for specific configurations.

In this study, we propose a hybrid numerical formulation for curvilinear grids based on Ducros's sensor [START_REF] Ducros | Large-eddy simulation of the shock/turbulence interaction[END_REF] that is able to switch from a low-dissipation formulation in presence of vortical structures to a robust formulation in presence of high gradients.

Organization of the paper

The paper is organized as follows: in Section 2 we describe the physical models chosen to conduct the study, their numerical implementation is then detailed in Section 3. The test case and the meshing strategy are presented in Section 4 and the computational results are finally analyzed in Section 5.

Physical modelling

Governing equations

To describe the compressible gas flows of interest in this study, one can use the general conservative form of the averaged/filtered Navier Stokes equations:

∂w c ∂t + div(F) = div(G) (1)
w c = t (ρ, ρV, ρE) denotes the vector of conservative quantities. ρ is the gas density, V = t (u, v, w) is the velocity vector of the flow and E represents the total energy of the gas. F and G denote respectively the convection and diffusion terms:

F = t (ρV, ρV ⊗ V + P I, (ρE + P )V) (2)
G = t (0, τ * , τ * .V -q * ) (3)
with P the thermodynamic pressure, τ * the combination of the viscous stress and the modelled turbulent stress (Reynolds or subgrid) tensors:

τ * ij = 2(µ + µ t )(S ij - 1 3 δ ij S ii ) i, j = 1..3 (4)
S ij = 1 2 ( ∂u i ∂x j + ∂u j ∂x i ) i, j = 1..3 (5)
with µ and µ t being the dynamic viscosity and the modelled turbulent viscosity which is provided by a turbulence model. q * = q + q t is the total heat flux with q the laminar heat flux computed with Fourier's law and q t the modelled turbulent heat flux linked to µ t by a turbulent Prandtl number.

Turbulence modelling with Zonal Detached Eddy Simulation (ZDES)

ZDES is a multiresolution approach developed by ONERA [START_REF] Deck | Recent improvements in the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF][43] which aims at treating in a simple model all classes of flow problems indicated in Fig 2 . In this formulation, three specific hybrid length scale formulations (see Eq. ( 6) ) also called modes are optimized to be employed on three typical flowfield topologies as illustrated in Fig

2 .
In the frame of the SA model, dZDES replaces the distance to the wall d w according to:

dZDES = dmode ∆, d w , U i,j , ν, ν, ∂ ν ∂n , ∂||ω|| ∂n (6) mode = 1, 2, 3 (7)
with ∆ denoting the subgrid length scale entering dZDES that depends on the grid cell size, U i,j the velocity gradients, d w the distance to the wall, ||ω|| the magnitude of vorticity, ν and ν the kinematic and pseudo eddy viscosity as well as ∂ ν ∂n the derivative in the wall-normal direction. The definition of dZDES is not a minor adjustment in the DES framework since the modified lengthscales depend not only on the grid but also on the velocity and eddy viscosity fields and their wall-normal derivatives. ZDES takes full advantage of its zonal nature since complex configurations can be accurately treated by the use of various formulations within the same calculation (an example is provided in [START_REF] Deck | Numerical investigation of the flow dynamics past a three-element aerofoil[END_REF]). To study generic axisymmetric base flows for which the separation is clearly fixed by the geometry, the mode 1 of the ZDES is an appropriate choice, it will be used in this study. However, when studying more realistic space launcher geometries as in [3][45], determining the separation point a priori can be difficult and the use of a completely automatic RANS/LES formulation such as the mode 2 of the ZDES is needed. To this end, the latest developement of the automatic mode of ZDES, refered to as ZDES mode 2 (2020) [START_REF] Deck | Towards an enhanced protection of attached boundary layers in hybrid RANS/LES methods[END_REF] will be assessed in this study as well.

Numerical Methods

Discrete formulation

We used the FLU3M solver [START_REF] Guillen | Design of a 3D multidomain Euler code[END_REF] developed by ONERA to solve the compressible Navier-Stokes equations on multiblock structured grids. An implicit finite volume form of the governing equations (Eq.( 1)), based on Gear's second-order accurate time discretization, yields for a hexahedral cell centered in (i,j,k) with a volume Ω ijk :

3 2 ∆ (wc) n+1 ijk - 1 2 ∆ (wc) n ijk + ∆t Ω ijk 6 l=1 (F(wc) -G(wc)) n ijk,l = - ∆t Ω ijk 6 l=1 ∂F(wc) ∂wc - ∂G(wc) ∂wc n ijk,l ∆(wc) n+1 ijk ( 8
)
Where n and n+1 represent consecutive time steps and ∆ (w c )

n+1 ijk = (w c ) n+1 ijk -(w c ) n ijk .
For each interface l associated to a normal vector n= t (n x , n y , n z ) and a surface ∆S ijk,l , the discrete convective and diffusive fluxes F and G are defined as:

F(wc) ijk,l =            ρ(V.n) ρu(V.n) + P nx ρv(V.n) + P ny ρw(V.n) + P nz (ρE + P )(V.n)            ijk,l .∆S ijk,l (9)
G(wc) ijk,l =            0 τ * xx nx + τ * xy ny + τ * xz nz τ * xy nx + τ * yy ny + τ * yz nz τ * xz nx + τ * yz ny + τ * zz nz τ * .V -q * .n            ijk,l .∆S ijk,l (10)
In order to use Eq.( 8) to compute (w c ) n+1 ijk at the center of cells grid knowing (w c ) n ijk , we thus have to choose numerical procedures to evaluate F(w c ) ijk,l , G(w c ) ijk,l and there Jacobian matrices ∂F(wc) ∂wc n ijk,l , ∂G(wc) ∂wc n ijk,l at the interface between cells.

To assess the Jacobian matrices we use the procedure detailed by Pechier in [START_REF] Pechier | Prévisions numériques de l'effet Magnus pour des configurations de munitions[END_REF] [START_REF] Pechier | Magnus effect over finned projectiles[END_REF]. Furthermore, to evaluate the diffusion flux vector G, values of ρ, D * , µ * , λ * ,V, h i at the cell interface are needed, the average of the right and left center values are used to estimate them. Estimations of velocity and temperature gradients are needed as well. To this end, a corrected Green-Gauss approach described in [START_REF] Manueco | On the coupling of wall-model immersed boundary conditions and curvilinear body-fitted grids for the simulation of complex geometries[END_REF] is used. It is designed to prevent numerical instabilities that can arise due to odd-even decoupling oscillations allowed by the original Green-Gauss approach.

Finally, two numerical features will be used to evaluate convective fluxes, namely a numerical scheme and a numerical procedure designed to increase the scheme accuracy: the MUSCL reconstruction method. On the basis of the literature review displayed in Sec.1.3, we propose here a hybrid numerical scheme formulation blending the low-dissipation AUSM+(P) [START_REF] Mary | Large Eddy Simulation of Flow Around an Airfoil Near Stall[END_REF] and the original AUSM+ [START_REF] Liou | A sequel to AUSM: AUSM+[END_REF] combined with different MUSCL reconstruction strategies. The hybridisation is based on the use of Ducros's sensor [START_REF] Ducros | Large-eddy simulation of the shock/turbulence interaction[END_REF] to detect zones featuring shock waves in order to add numerical dissipation locally. The idea here is to preserve the good scale resolution properties of the low-dissipation AUSM+(P) scheme in vortical regions and to enhance its robustness near shocks. The different aspects of this framework are detailed below.

Numerical scheme

Low-dissipation AUSM+(P)

Mary and Sagaut [START_REF] Mary | Large Eddy Simulation of Flow Around an Airfoil Near Stall[END_REF] proposed a simplified formulation of the AUSM+(P) designed to ensure a low numerical dissipation. We recall here the main features of this scheme. For the sake of explanation, we suppose that we are willing to compute this low dissipation (LD) flux F LD j+ 1 2 through the interface between the cell (i,j,k) and the cell (i,j+1,k). The Euler fluxes are approximated by :

F LD j+ 1 2 = V 1 Ψ L + Ψ R 2 -|U dis | Ψ R -Ψ L 2 + P 1 (11)
With Ψ L,R being the left and right states obtained at the interface with the chosen MUSCL reconstruction (described in Sec.3.3). V 1 corresponds to the interface fluid velocity and

P 1 = [0, P L +P R 2 n x , P L +P R 2 n y , P L +P R 2 n z , 0]. U dis is a
parameter based on the local fluid velocity to control the numerical dissipation. More precisely, V 1 is computed with a pressure stabilization term:

V 1 = (V L + V R ).n 2 -c 2 (P R -P L ), (12)
and U dis is defined as: ζ ∈ {0; 1} is a sensor used to minimize dissipation by enabling a fully centered scheme in monovariant zones [START_REF] Mary | Large Eddy Simulation of Flow Around an Airfoil Near Stall[END_REF].

U dis = ζ * max(| (V L + V R ).n 2 |, c 1 ) (13)
As explained earlier, this scheme has been used with success on various cases but is not adapted for flow configurations leading to strong shocks. For the case with propulsive jet studied here, this scheme, used with a 3rd-order MUSCL reconstruction, leads to a rapid numerical divergence. This is also the case when:

The scheme is used with a 2nd-order MUSCL reconstruction.

The scheme is used with ζ forced to 1 and with a 2nd-order MUSCL reconstruction.

As explained in Sec1.3, we want to enhance locally the dissipation of this scheme. One idea is to blend this formulation with another numerical scheme providing good shock capturing properties. To this end we chose the AUSM+ formulation described here after.

AUSM+

The AUSM+ scheme has been proposed by Liou [START_REF] Liou | A sequel to AUSM: AUSM+[END_REF]. The main idea here is to decompose the non-viscous fluxes F(w c ) between convective fluxes F (c) and acoustic fluxes P.

F(w c ) = F (c) + P (14)
F (c) = M a t (ρ, ρV, ρH), P = t (0, P x , P y , P z , 0) (15)
M= V.n a , a the speed of sound and H the total enthalpy (H = E + P ρ ). We need to express a convective flux f (c) j+1/2 and a pressure flux p j+1/2 . Defining a j+1/2 as the speed of sound at the cell interface, one can express:

f (c) j+1/2 = m j+1/2 a j+1/2 Ψ j+1/2 , Ψ = t (ρ, ρV, ρH) (16)
and

p j+1/2 = t (0, P x , P y , P z , 0) j+1/2 (17)
Ψ j+1/2 is given by a basic upwinding procedure:

Ψ j+1/2 =    Ψ L if m j+1/2 ≥ 0 Ψ R otherwise (18)
The terms m j+1/2 , p j+1/2 and a j+1/2 have yet to be defined. Liou [START_REF] Liou | A sequel to AUSM: AUSM+[END_REF] introduces a decomposition of these terms as the sum of two contributions:

m j+1/2 = M + (M j ) + M -(M j+1 ) (19)
p j+1/2 = P + (M j )p j + P -(M j+1 )p j+1 (20)
Liou [START_REF] Liou | A sequel to AUSM: AUSM+[END_REF] proposes the forms:

M ± (M ) =      1 2 (M ± |M |), if |M | ≥ 1 ± 1 2 (M ± 1) 2 ± 1 8 (M 2 -1) 2 otherwise (21)
P ± (M ) =      1 2 (1 ± sign(M )), if |M | ≥ 1 ± 1 4 (M ± 1) 2 (2 ± M ) 2 ± 3 16 M (M 2 -1) 2 otherwise (22)
with sign(M ) = M |M | . For the speed of sound at the interface, a simple average yields satisfying results:

a j+1/2 = 1 2 (a j + a j+1 ) (23)
The AUSM+ numerical flux F AU j+ 1 2 is finally computed as:

F AU j+ 1 2 = m j+1/2 a j+1/2 Ψ j+1/2 + p j+1/2 (24)

Hybrid shock capturing scheme

To ensure both a low numerical dissipation and good local shock capturing properties, we propose the following hybrid formulation:

F j+ 1 2 = (1 -β)F LD j+ 1 2 + βF AU j+ 1 2 (25)
With β a sensor designed to detect shocks, based on Ducros's [START_REF] Ducros | Large-eddy simulation of the shock/turbulence interaction[END_REF] sensor Φ:

Φ = (∇.u) 2 (∇.u) 2 + (∇ ∧ u) 2 + = 1.10 -32 s -2 (26)
∇.u is the divergence of the velocity vector, it takes high values near shocks and ∇ ∧ u evaluates the vorticity of the flow, it takes high values in vortical regions. Thus, Φ ∈ [0; 1] and it takes values close to one near a shock wave and close to zero in vortical regions. A more complete formulation has been derived to enhance dissipation near shocks, to maintain β between 0 and 1 and to ensure a smooth variation of β limiting sensitivity to numerical instabilities:

β(Φ) = tanh(c * Φ) (27)
As shown in Fig. 3, the choice of c=3 for the constant value provides a good compromise. Indeed, this value ensures a high β value plateau near shock waves (Φ close to 1) while maintaining low values of β(Φ) in vortical region where usually Φ ≤ 0.01, for instance β(0.01) ≈ 0.03.

MUSCL reconstruction strategy 185

To complete the description of the numerical scheme presented above, the procedure providing the left and right states Ψ L,R at the interface has to be clarified.

Slope computation

As explained in [START_REF] Koren | A robust upwind discretization method for advection, diffusion and source terms[END_REF][52] [START_REF] Chauvet | Simulation numérique et analyse physique d'un jet propulsif contrôlé par des injections radiales[END_REF], if one simply uses the center values Ψ L = (w c ) n i,j,k and Ψ R = (w c ) n i,j+1,k to compute F j+ 1 2 , the presented numerical scheme would only be first order accurate in space. In order to enhance the scheme accuracy, Van Leer [START_REF] Van Leer | MUSCL, a new approach to numerical gas dynamics[END_REF] has proposed the MUSCL strategy (Monotonic Upstream-centered Scheme for Conservation Law) which consists in considering a linear variation of state variables inside cells instead of using constant values as shown in Fig. 4. For a variable w of the vector (w c ), two types of evaluations can be used to compute the slope of variations inside the cells:

U pwind :      w L j+ 1 2 = w j + w j -w j-1 2 w R j+ 1 2 = w j+1 - w j+2 -w j+1 2 (28)
Centered :

     w L j+ 1 2 = w j + w j+1 -w j 2 w R j+ 1 2 = w j+1 - w j+1 -w j 2 (29)
The general framework of the MUSCL reconstruction method consists in a combination of the two mentioned methods thanks to a parameter κ : with

     w L j+ 1 2 = w j + 1 -κ 4 ∆ j-1 2 + 1 + κ 4 ∆ j+ 1 2 w R j+ 1 2 = w j+1 - 1 + κ 4 ∆ j+ 1 2 - 1 -κ 4 ∆ j+ 3 2 (30)
∆ j+ 1 2 = w j+1 -w j ; ∆ j-1 2 = w j -w j-1 ; ∆ j+ 3 2 = w j+2 -w j+1 (31)
The upwind formulation is obtained with κ = -1, a centered formulation (2nd-Order) is obtained with κ = 1 and a third-order upwind biased (denoted as UB) formulation is obtained with κ = -1 3 [START_REF] Chauvet | Simulation numérique et analyse physique d'un jet propulsif contrôlé par des injections radiales[END_REF].

centered formulation (κ = 1) combined with the Minmod limiter [START_REF] Sweby | High resolution schemes using flux limiters for hyperbolic conservation laws[END_REF] Γ :

     w L j+ 1 2 = w j + 1 2 Γ(∆ j+ 1 2 , ∆ j-1 2 ) w R j+ 1 2 = w j+1 - 1 2 Γ(∆ j+ 1 2 , ∆ j-1 2 ) (32)
The Minmod limiter thus selects the smallest slope in monovariant zones and reduces to first order accuracy around local extrema. The less dissipative strategy is to use the κ = -1/3 slope reconstruction without limiter, namely the third-order Upwind Biased formulation (called UB). Finally, we also use a limited version of the UB formulation thanks to a Koren-type limiter K(r) presented in [START_REF] Collercandy | An improved approach for the computation of transonic/supersonic flows with applications to aerospace configurations[END_REF] :

200          w L j+ 1 2 = w j + K( ∆ j+ 1 2 ∆ j-1 2 ) 1 -κ 4 ∆ j-1 2 + 1 + κ 4 ∆ j+ 1 2 w R j+ 1 2 = w j+1 -K( ∆ j+ 3 2 ∆ j+ 1 2 ) 1 + κ 4 ∆ j+ 1 2 - 1 -κ 4 ∆ j+ 3 2 (34)
With κ = -1 3 and:

Koren : K(r) = 3r + 2 1 + r + r 2 + 2 , 2 = 1e -32 (35)
As showed in Fig. 5, this formulation (called Koren here) is close to the third-order Upwind Biased one (UB)

when K(r) ≈ 1, that is to say when the ratio r between the right and left slopes is around 1 and is limited otherwise

K(r) < 1.

Hybrid MUSCL formulation 205

As described in Sec. 1.3, the good numerical accuracy showed in [START_REF] Weiss | On the dynamics of axisymmetric turbulent separating/reattaching flows[END_REF] [17] is obtained thanks to the combination of a low dissipation AUSM+(P) scheme and a UB MUSCL reconstruction. The hybrid scheme formulation presented above enables to retain the accuracy of the low dissipation AUSM+(P) in desired regions and the selection of the MUSCL reconstruction method is discussed in the following. Indeed, it is shown in [START_REF] Reynaud | Towards binary gas mixture ZDES for space launcher base flow prediction[END_REF] that the AUSM+ scheme, which is used for the shock capturing part of the hybrid scheme, enables to compute strong shock tube cases with a Minmod or a Koren reconstruction but not with the UB reconstruction. We thus propose here a hybrid MUSCL reconstruction strategy to complete the hybrid numerical scheme formulation. At a cell interface the left and right states are computed as:

Ψ L,R = (1 -β )Ψ (3)
L,R + β Ψ (2) L,R (36)
With Ψ

L,R being the state obtained with a 2nd-order MUSCL reconstruction (here Minmod) and Ψ

L,R being the state obtained with a 3rd-order MUSCL reconstruction (here Koren or UB). β is equal to β (see Eq.27) but can be forced to 1 in desired zones.

Based on this framework we have used two different reconstruction strategies for computations with jet. First, a Koren/Minmod strategy is tested with a Koren reconstruction for Ψ L,R , β = β out of the jet zones and β = 1 in the jet zones. Indeed if β is not forced to 1 in jet zones, that is to say if a full Minmod reconstruction is not used, the UB/Minmod reconstruction is not robust enough and the computation diverges.

Test case

The test case configuration is an axisymmetric afterbody studied experimentally by Depres et al. ( [START_REF] Deprés | Analysis of unsteadiness in afterbody transonic flows[END_REF], [START_REF] Deprés | Analyse physique et modélisation des instationnarités dans les écoulements d'arrière-corps transsoniques[END_REF]) in ONERA's S3Ch wind tunnel. Experiments were performed with and without a cold air propulsive jet.

Experimental setup

As shown in Fig. 6a The free stream Mach number is set to 0.702 in the wind tunnel and the nozzle pressure ratio (NPR= Pt j P0 ) of total jet pressure (P tj ) over freestream pressure (P 0 ) is 34.2, which corresponds to an adapted jet i.e.(exit jet pressure close to the external pressure). The Reynolds number of the external flow, based on the forebody diameter, is Re D is 1.2 × 10 6 .

Meshing

The mesh is presented in Fig. (6b). In coherence with experimental measurements, the computed forebody length was chosen to ensure a ratio δ D = 0.2 between the external boundary layer thickness δ and the forebody diameter. 40 points are clustered in this boundary layer and the ∆y + = 1 condition is satisfied. In Fig. 7, one can verify that the thickness of the computed boundary layer compares well with experimental measurements for different azimuthal angles. The multi-block structured meshes used contain 12.5M cells for the no-jet case and 16.2M cells for the case with jet (Fig. 6c). 240 points are distributed in the azimuthal direction resulting in a discretization ∆θ = 1.5 • . Besides, an O-H grid topology is adopted to avoid a singularity on the axes.

Numerical setup

Ten different computations have been performed for this study to compare the several methodologies described in the previous sections. For each of the configurations (with and without jet), four different ZDES computations and a RANS AUSM+ computation have been run, as summarized in table 3. ZDES computations were run with 5 sub-iterations. First, we let the flow develop for 200 ms (≈ 100 × T s with Ts the vortex shedding period) to avoid The various computational configurations listed in Tab.3 enable a gradual assessment from the current numerical standards for a generic transonic afterbody geometry without jet (NJ1) to our proposition for a numerical framework adapted to complex space launcher geometries with propulsive jet (J4). During this process, comparing NJ1 vs NJ3 enables to evaluate the performance of the hybrid scheme (used in NJ3) compared to the LD AUSM+(P) scheme (used in NJ1) on a configuration without shocks. Then, the performance of the hybrid scheme for cases

Results

We present here the results obtained for the afterbody computations with and without jet. In addition to the physical description of the flow, this section will provide an evaluation of the developed numerical workflow. Indeed, numerical data are compared with both experimental measurements and previous numerical studies following the aforementioned nomenclature for validation of simulation techniques [START_REF] Weiss | On the coupling of a zonal body-fitted/immersed boundary method with ZDES: Application to the interactions on a realistic space launcher afterbody flow[END_REF] (see Tab 1).

Flow topology

Instantaneous visualisations

The salient features of the instantaneous flow described in [START_REF] Weiss | On the dynamics of axisymmetric turbulent separating/reattaching flows[END_REF] are evidenced in Fig. 8 and Fig. 9 by plotting an iso-surface of the Q-criterion. The roll-up of the mixing layer induces the formation of azimuthal coherent structures which are disrupted by azimuthal instabilities and progressively bent in the flow direction to form hairpin vortices as they are convected near the reattachment point. Passed the reattachment we mainly observe longitudinal eddies.

Note that almost no differences can be seen between NJ1 and NJ4 computations showing that the desired properties of the proposed hybrid scheme in vortical regions are achieved. A sliced view of the flow provided in Fig. 10 confirms the global similarity between the two computations but enables to observe that the NJ1 computation resolves more turbulent structures close to the base. In other words, the use of the hybrid scheme and the Koren limiter induces a minor increase of the numerical dissipation in smooth regions. The presence of the jet tends to amplify the stretching and the alignment of the structures with the flow direction after the extension due to a higher entrainment effect.

Furthermore, for the case with jet, the J3 computation (Koren/Minmod) allows to advect more structures from the recirculation area into the jet area than the J1 computation (Minmod in the jet area). This observation is confirmed in Fig. 10. The J3/J4 hybrid MUSCL strategy is thus both less dissipative in the jet area and easier to use than the J1/J2 MUSCL strategy. Indeed, the decision load of the user is very low for computation J4 which is automatic on both numerical and turbulence treatment aspects. Fig. 10 also displays the strong normal shock structure of the jet which is classically the source of numerical difficulties. Overall, the present ZDES computations using the hybrid scheme are robust and permit a noticeable improvement in scale resolving compared to results obtained in [START_REF] Deck | Detached and large eddy simulation of unsteady side-loads over an axisymmetric afterbody[END_REF] using Roe's scheme on a similar mesh.

Mean flow topology

The main properties of the mean velocity fields are shown in Fig. 11. Let us be reminded that for a short extension length, the flow does not reattach on the solid part of the rear body (see Fig. 1). No PIV measurements were made during Depres' experiments [START_REF] Deprés | Analysis of unsteadiness in afterbody transonic flows[END_REF] but PIV results from two other experimental campaigns [START_REF] Van Gent | From particles to pressure PIV-based pressure reconstruction for base flows[END_REF][24] can be used here for the analysis. The test conditions during these campaigns, [START_REF] Van Gent | From particles to pressure PIV-based pressure reconstruction for base flows[END_REF] (D=50 mm, d/D=0.34, L/D=0.6, M ∞ = 0.76, model held by a strut, tests with and without jet) and [START_REF] Saile | Experimental analysis on near-wake flows of space transportation systems[END_REF] (D=66.7 mm, d/D=0.4, L/D=0.6, M ∞ = 0.8, tests with jet), were close but nevertheless not equivalent to the present cases.

We observe in Fig. 12 that the presence of the jet causes both a shortening of the recirculation bubble and an increase of the backflow velocity around the extension. This effect is predicted by both RANS and ZDES computations.

More precisely, in the absence of the jet, ZDES computations predict a longer recirculation bubble ( Lr D ≈ 1.5) than the RANS computation ( Lr D ≈ 1.16). For this particular flow characteristic, the RANS predictions appear closer

Statistical properties of the wall pressure field

The distribution of mean (Cp) and fluctuating (Cp rms ) pressure coefficients on the base and on the extension are discussed here. Fig. 14 indicates an important decrease in base pressure (i.e. an increase in base drag) due the presence of the jet as shown by both experimental and numerical results. Indeed, ZDES calculations for the jet-off case exhibit a flat Cp profile with -0.138 ≤ Cp ≤ -0.135 illustrating the good agreement with experimental results 300 in which -0.131 ≤ Cp ≤ -0.124. However, when jet is on, the Cp levels fall to -0.18 ≥ Cp ≥ -0.185 according to ZDES computations and -0.195 ≥ Cp ≥ -0.21 on wind tunnel measurements. The observed discrepancy in this case points out that the aspiration effect of the jet is underestimated numerically. RANS results present major discrepancies both in mean pressure levels and on their radial variation.

The pressure drop induced by the jet is even more important when looking at the longitudinal variation of Cp on the extension presented in Fig. 15. One can notice a global improvement on pressure level predictions with current ZDES computations compared to [START_REF] Deck | Detached and large eddy simulation of unsteady side-loads over an axisymmetric afterbody[END_REF] both with and without the jet. Present RANS results are very far from being in agreement with the experimental data. This is consistent with the overestimation of the recirculation velocity previously described. Similar disparities for classical RANS models have been reported in [START_REF] Simon | Zonal-Detached-Eddy Simulation of Projectiles in the subsonic and transonic regimes[END_REF] [START_REF] Catalano | Cfd contribution to the aerodynamic data set of the vega launcher[END_REF].

Focusing on the distribution of pressure fluctuations along the extension (Fig. 16), we observe a continuous increase in the intensity of the fluctuations in the x-direction both with and without the jet. Furthermore, as described in the literature, the overall Cp rms levels are noticeably higher in the presence of the jet, especially near the end of the extension. The fluctuation levels are slightly underestimated by ZDES computations near the base but their rise in the second half of the extension is well captured. Let us be reminded that the incoming attached boundary layer is treated in URANS within mode 1 and mode 2 of ZDES. The assessment of incoming fluctuations inherited from the boundary layer is beyond the scope of this study but could be envisaged in the future within mode 3 of ZDES [START_REF] Deck | Zonal Detached Eddy Simulation (ZDES) of a spatially developing 575 flat plate turbulent boundary layer over the reynolds number range 3 150 ≤ Re θ ≤ 14 000[END_REF].

For mean values and second order statistics of wall pressure, the displayed results show that the different ZDES computations provide results very close to each other. For the case without jet, a minor discrepancy is the behaviour of the mode 2 computation (NJ2) at the end of the extension. Indeed, compared to the other mode 1 results, it slightly underestimates the Cp rms levels (Fig. 16) and overestimates the mean pressure level (Fig. 15) for x/D ≥0.45.

This difference can be attributed to the deactivation of wall functions (used in RANS and mode 2 of the ZDES) when using the mode 1 of the ZDES (see [START_REF] Deck | Recent improvements in the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF]), which induces a different treatment of the small reattachment bubble observed at the end of the extension on Fig. 12. In the presence of the jet, one can observe that the two computations using the Koren/Minmod reconstruction (J3 and J4) tends to slightly underestimate the level of 325 fluctuations (Fig. 16) in comparison with the UB/Minmod computations (J1, J2). This can be attributed to the extra dissipation induced by the use of the Koren limiter compared to the unlimited Upwind Biased slope in the base area as seen in Fig. 10. However, the two computations (J3 and J4) using a hybrid Ψ (3) /Ψ (2) limiter formulation (see Sect. 3.3.3) in the jet area found a lower pressure coefficient near the jet compared to the two using a Ψ (2) limiter in this area (J1; J2). In this zone, the formulation used in (J3; J4) is the less dissipative and provides an 330 increase of the jet aspiration effect. This effect on pressure levels is nevertheless limited to the area near the jet and does not provide a decrease on base pressure levels.

Spectral analysis

To study more deeply the dynamic behaviour of the base flow, one can examine the spectral distribution of the pressure fluctuations mentioned above to determine characteristic frequencies. Indeed, Cp rms = prms q∞ with q 2 ∞ = 1 2 γp ∞ M 2 ∞ the dynamic pressure and one can express p rms values as: shown in the broadband spectra in Fig. 18. Such a contribution is found numerically, even though it is underestimated compared to the experiment. Indeed, as regards base pressure fluctuations, the level of fluctuations is underestimated by all ZDES computations. The smaller duration of the pressure signal and the absence of turbulent fluctuations within the incoming boundary layer can again partially explain those discrepancies but it is also worth noticing that the experimental probes (r/R=0.7) are located just under the delimitation of the recirculation bubble found numerically (Fig. 13). The discrepancies observed in pressure fluctuations on the base could thus be partially caused by a different relative position between the probe and the recirculation bubble in the experiment.

p rms = ∞ 0 G(f )d(f ) = ∞ 0 G(St D )d(St D ) (37)
Apart from the slight underestimations of the overall fluctuation level found in Fig. 16 with the J3 and J4 computations and again noticeable in Fig. 17 and Fig. 18 , the different ZDES computations again display very close results. This indicates that the automatic shock-capturing formulation proposed (J4) is able to retain most of the shock-free reference framework (NJ1) accuracy. This is further confirmed by the two-point spectral analysis displayed below.

Two-point analysis

To study the spatial organization of the flow at the previously highlighted frequencies, the azimuthal coherence of two pressure sensors p 1 (r, x, φ 1 ) and p 2 (r, x, φ 2 ) located in a plane normal to the inflow at a constant position

x/D and a constant radius r/R can be considered. Assuming the hypothesis of an homogeneous flow, i.e. without any preferred angle of reference φ 1 , the complex coherence function may be expressed as:

C(f, r, x, ∆φ) = (C r + jC i )(f, r, x, ∆φ) = S 12 (f, r, x, ∆φ) S 1 (f, r, x, ∆φ)S 2 (f, r, x, ∆φ) (38)
where j = √ -1, and C r and C i are the real and imaginary part of the cross-spectral density function S 12 and ∆φ = φ 1 -φ 2 . Assuming that the disturbances do not exhibit any particular direction of propagation, one has S 12 (∆φ) = S 12 (-∆φ). Besides, the hypothesis of isotropy yields C i = 0. Consequently, the C r function is 2πperiodic with respect to ∆φ and can, thanks to a Fourier transform in azimuthal modes, be expressed as:

C r (f, ∆φ) = ∞ m=0 C r,m (f )cos(m∆φ). ( 39

Conclusion

Starting from a standard, non-automatic, numerical framework for accurate RANS/LES computations of shockfree space launcher afterbody flows, we have gradually built an efficient automatic workflow adapted to turbulent flows containing shocks. This strategy is based on several numerical features. First, the implementation of a hybrid numerical scheme and a hybrid MUSCL reconstruction based on Ducros's sensor [START_REF] Ducros | Large-eddy simulation of the shock/turbulence interaction[END_REF] enables to switch from a lowdissipation formulation in presence of vortical structures to a robust formulation in presence of high gradients. Then, the latest automatic mode of ZDES (ZDES mode 2 (2020) [START_REF] Deck | Towards an enhanced protection of attached boundary layers in hybrid RANS/LES methods[END_REF]) is used for turbulence modelling. The contributions of these features have been progressively evaluated thanks to eight ZDES computations of a transonic space launcher afterbody geometry with and without propulsive jet. Thorough comparisons with experimental measurements and standard ZDES computations have shown that, notwithstanding a minor increase of numerical diffusion in vortical regions, the proposed automatic shock-capturing RANS/LES workflow is able to retain the most of the shock-free reference framework accuracy and provides satisfying flow predictions according to an extended nomenclature for validation of simulation techniques [START_REF] Weiss | On the coupling of a zonal body-fitted/immersed boundary method with ZDES: Application to the interactions on a realistic space launcher afterbody flow[END_REF]. The presented workflow thus appears well suited to accurately study realistic launcher afterbody geometries with propulsive jet(s) with a reasonable computational cost and a low user workload.

As regards the overall results provided by the ZDES computations, the introduction of turbulent content within the space launcher main body boundary layer (e.g. mode 3 of the ZDES [START_REF] Deck | Zonal Detached Eddy Simulation (ZDES) of a spatially developing 575 flat plate turbulent boundary layer over the reynolds number range 3 150 ≤ Re θ ≤ 14 000[END_REF]) and the taking into account of the thermodynamic properties of the propulsive jet could be the next steps to further improve the prediction of pressure levels for realistic space launchers.

Figure 1 :

 1 Figure 1: Topology of the interactions between a propulsive jet and a subsonic base flow

Figure 2 :

 2 Figure 2: Classification of typical flow problems. I: separation fixed by the geometry, II: separation induced by a pressure gradient on a curved surface, III: separation strongly influenced by the dynamics of the incoming boundary layer

c 1 and c 2

 2 are constant parameters based on the flow reference values (c 1 = 0.04 * U ∞ and c 2 = 0.04 * 1 U∞ρ∞ ) and

Figure 3 :

 3 Figure 3: Evolution of function tanh(c * Φ) for different values of c.

Figure 4 :

 4 Figure 4: Schemes of the upwind (left) and centered (right) slope evaluation method.

  with M inmod : Γ(a, b) = sign(a)max(0, min(|a|, sign(a)b))

Figure 5 :

 5 Figure 5: Variation of the Koren limiter.

  and β = β. Second a UB/Minmod strategy is assessed with a 3rd-Order upwind biased (UB) reconstruction for Ψ

  the model consists in an axisymmetric blunt body with a diameter D=2R equal to 100 mm and a 60mm-long extension with a smaller diameter d=2r equal to 40 mm. The ratios L D = 0.6 and d D = 0.4 are representative of current space launcher designs [57] [5].

Figure 6 :

 6 Figure 6: Test case setup

  240

Figure 7 :

 7 Figure 7: Boundary layer profile on the forebody at the station X/D=-2.45

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: Instantaneous flow visualization of the NJ1 (left) and NJ4 (right) computations for the case without jet. Iso-surfaces of Q criterion (Q*D 2 /U 2 ref =50) and iso-contours of pressure coefficient. (See Tab. 3 for computation details)

Figure 11 :Figure 12 :

 1112 Figure 11: Mean flow topology for RANS and ZDES NJ1 for the case without jet (left) and RANS and NJ1 for the case with jet (right). isoline u U ref =0. (See Tab. 3 for computation details)

Figure 13 :

 13 Figure 13: Mean flow topology with the four ZDES set-ups for the case with jet.isoline u U ref =0. (See Tab. 3 for computation

Figure 14 :

 14 Figure 14: Distribution of the pressure coefficient on the base without jet (top) and with jet (bottom). (See Tab. 3 for computation details)

Figure 15 :

 15 Figure 15: Distribution of the pressure coefficient on the extension without jet (top) and with jet (bottom). (See Tab. 3 for computation details)

Figure 16 : 2 ∞.Fig. 17

 16217 Figure 16: Distribution of the pressure fluctuations on the extension without jet (top) and with jet (bottom). (See Tab. 3 for computation details)

Fig. 18 shows

 18 Fig.18shows the pressure power spectral density on the base of the afterbody. If the presence of a peak around St D ≈ 0.2 is found both experimentally and numerically for the case without jet, showing that the vortex shedding phenomenon also affects the pressure fluctuations on the base, no clear sign of this contribution is found in the presence of the jet. In both cases, the main source of pressure fluctuations is a low frequency range St D ≤ 0.1 as

  ) C r,m represents the percentage of the fluctuating energy at frequency f relative to the azimuthal mode m since C r,m = 1. Let us be reminded that m = 0 and m = 1 modes are characterized by an inphase and antiphase

Figure 17 :Figure 18 :

 1718 Figure 17: Power spectral density of the fluctuating pressure near the end of the extension (x/D=0.55) without jet (top) and with jet (bottom). (See Tab. 3 for computation details)

Figure 19 :Figure 20 :

 1920 Figure 19: Spectra of the first 3 azimuthal modes for the case without jet. (a)-(c)-(e) on the extension (x/D=0.55). (b)-(d)-(f) on the base (r/R=0.7). (See Tab. 3 for computation details)

Table 2 :

 2 Overview of studies on high subsonic base flows over axisymmetric generic geometries. Exp : Wind-tunnel experiments.

	Adapted from [22]

Table 3 :

 3

List of computational setups used. UB= Upwind Biased (κ = -1 3 in Eq.30), LD= Low Dissipative version

[START_REF] Mary | Large Eddy Simulation of Flow Around an Airfoil Near Stall[END_REF]

.

of the extension. This observation can be related to PIV results[START_REF] Van Gent | Effects of exhaust plume and nozzle length on compressible base flows[END_REF] showing that the probability of reversed flow

Acknowledgments

The authors wish to thank the Centre National d' Études Spatiales (CNES) for financial support. The Ph.D. work of J. Reynaud is funded by CNES and ONERA.

Slope limiters

The extension of the spatial accuracy of the numerical scheme can lead to the apparition of spurious oscillations in the presence of high gradients. To avoid these numerical instabilities degrading the robustness of the computation, the slopes presented in the previous section have to be limited. In this study three MUSCL strategies are used and mixed as described in the next section. The most dissipative strategy (called Minmod) consists in the use of a