domain: spi phys

text:
Vincent Drouet
email: vincent.drouet@onera.fr

Ysolde Prévereaud
email: ysolde.prevereaud@onera.fr

Jean-Marc Moschetta
email: jean-marc.moschetta@isae-supaero.fr

Nathalie Bartoli
email: nathalie.bartoli@onera.fr

Sylvain Dubreuil
email: sylvain.dubreuil@onera.fr

Julien Annaloro
email: julien.annaloro@cnes.fr

Reduced order models for heat flux and pressure distributions on space debris afterbodies

Keywords: orbital debris, aerothermodynamics, hypersonic flows, afterbodies, modeling

The prediction of the on-ground risk caused by re-entering space debris requires accurate and computationally affordable models to compute wall heat flux and pressure coefficient C p at every stage of the reentry. Such models already exist for the wind area of the debris, i.e. the area directly impinged upon by fictitious lines parallel to the incoming flow. But the heat flux and C p are often neglected in the shadow area, even though they can be relatively high for certain shapes of debris. The present work focuses on developing reduced order models for heat flux and C p distributions in the shadow area of space debris, for hypersonic continuous flow conditions. We identified four phenomena which appear in the shadow area and cause relatively high levels of these aerodynamic quantities (attached flow, detached flow with fluid reattachment, detached flow with solid reattachment and shock-shock interactions). Models were developed for the heat flux and C p distributions caused by attached flows on the lee side of cylindrical geometries using the Proper Orthogonal Decomposition (POD) and interpolation method. Using this method, it is possible to reduce the number of required data by efficiently exploring the parameters domain of variation. The sample points were chosen thanks to an adaptive design of experiments, and the input data for the models were obtained by 3D Navier-Stokes computations of the flow around cylinders at incidence. The analysis of the computational results highlighted the influence of the Reynolds number based on the cylinder diameter Re D,∞ on the heat flux and pressure distributions. Reduced-order models were created for three input parameters (length L, diameter D and angle of attack α) and fixed incoming flow conditions, using the POD and interpolation method, and then were extended to other upstream hypersonic continuous flow conditions by more classical approaches based on non-dimensional quantities. Finally, the reentry trajectories of two different cylindrical debris were computed with and without the new models developed, to demonstrate their influence

Introduction

The uncontrolled atmospheric reentry of the Tiangong 1 Chinese space station in April 2018 highlighted the concern of on-ground risk posed by orbital debris. Currently, about 15 000 debris larger than 1 cm are in Earth orbit [START_REF] Anz-Meador | Orbital debris quarterly news[END_REF]. On average, one to two of these debris are subject to an uncontrolled atmospheric reentry each day, and one to two debris larger than 1 m are subject to an uncontrolled reentry each week. These numbers are likely to increase as several satellites constellations, formed of several hundreds or thousands of satellites, are being deployed [START_REF] Del Portillo | A technical comparison of three low earth orbit satellite constellation systems to provide global broadband[END_REF][START_REF] Pardini | Environmental sustainability of large satellite constellations in low earth orbit[END_REF]. The total casualty area forecast becomes a major issue for all space actors and especially for CNES, which is in charge of ensuring the strict application of the French Space Operation Law (LOS) by 2021, for French satellite-and launcher-operators and for launch operations from the French Guyana spaceport. Among others, this law requires that the maximum probability to have at least one victim does not exceed 1 × 10 -4 .

The on-ground risk posed by a reentering debris is determined by the geographical position and the size of the area of impact, and by the number and kinetic energy of the surviving fragments. These results are obtained by computing the trajectory, mass ablation and fragmentation of the debris, using engineering models for aerodynamic forces and heat flux on the walls of the debris at each timestep of the reentry. Depending on the modeling approach used, the software predicting the on-ground risk posed by orbital debris can be classified into two types, "spacecraft oriented" and "object-oriented" codes. The present study focuses on spacecraft oriented codes, like the software ARES developed at ONERA [START_REF] Prévereaud | Numerical and experimental study of the thermal degradation process during the atmospheric re-entry of a TiAl 6 V 4 tank[END_REF], which take into account the whole geometry of the debris before fragmentation and can compute the reentry trajectory with 6 degrees of freedom. Engineering models are used for the wall distributions of heat flux, pressure coefficient C p and friction coefficient C f on the whole surface of the debris. For the continuous hypersonic part of the reentry, the C f is negligible, the C p distribution is often predicted using Newton-like methods, and the wall heat flux is computed from C p using empirical relationships. These methods are only valid for the windward walls, i.e. the walls directly impacted by fictitious lines parallel to the freestream flow. In the "shadow region", defined by opposition to the windward walls and illustrated in figure 1, the C p is often set to 0, which results in a very low heat flux.

However, some phenomena described in section 2 can cause relatively high levels of C p and heat flux in the shadow region of space debris. For instance, studies have shown that the heat flux in the shadow region of a reentry body could reach 20 % of the heat flux at the stagnation point [START_REF] Wright | Aerothermal analysis of the project Fire II afterbody flow[END_REF][START_REF] Horvath | Blunt body near wake flow field at Mach 10[END_REF]. In such cases, neglecting heat flux distribution in the shadow region could lead to an overestimation of the risk posed by the debris, especially in the case of composite walls, for which the degradation process can occur at relatively low temperature. Furthermore, during the atmospheric reentry, debris can Figure 1: Definition of the "shadow region" (hatched zone) on a sphere; U ∞ is the direction of the freestream flow have a random tumbling motion. Thus, some walls of the debris can alternatively be in the windward region and in the shadow region. Neglecting the wall heat flux in the shadow region will result in a wrong prediction of the wall temperature, which can alter the prediction of the heat flux once the wall lies in the windward region again. Additionally, underestimating the C p levels in the shadow area can lead to overestimating the lift and drag coefficients. Navier-Stokes computations detailed in section 4 have shown that neglecting the C p distribution in the shadow area of a cylinder could lead to an error on the lift coefficient up to 180 %, which would result in a wrong prediction of the area of impact. Overestimating the drag coefficient also results in an underestimation of the debris speed, and therefore of its thermal degradation. Finally, a miscalculation of the C p distribution in the shadow region can lead to wrong values of the aerodynamic moments exerted on the debris, and to a wrong prediction of the debris orientation during the reentry, which can strongly affect the risk computation.

Therefore, the aim of this work is to provide more accurate models for wall heat flux and C p in the shadow region of space debris, during the hypersonic continuous phase of their atmospheric reentry. To do so, we identified the flow topologies causing relatively high levels of wall C p or heat flux in the shadow region of space debris, and present them briefly in section 2. The rest of the paper focuses on the modeling of two of these phenomena: the attached flows and the detached flows with fluid reattachment. In this work, we decided to model the wall heat flux and C p distributions on the lee side and on the aft of cylinders with different sizes (length L, diameter D and edge radius R e ), for different flight points (altitude h, Mach M ∞ and Reynolds numbers Re ∞ ) and angles of attack α. The variation range of input parameters is described in section 3. In section 4, we present the methodology used to build the models, and especially the Proper Orthogonal Decomposition (POD) and interpolation method, presented in section 4.1. This method requires an adapted design of experiments, presented in section 4.2. The entry data for the models was obtained by computational fluid dynamics (CFD) simulations, carried out with the reactive Navier-Stokes solver CEDRE from ONERA [START_REF] Refloch | CEDRE Software[END_REF][START_REF] Onera Cedre | [END_REF]. We present the physical models and the meshes used in the computations in section 4.3. The resulting wall C p and heat flux distributions in the shadow region of cylinders exhibited unpredicted behaviors and are detailed in section 5. Finally, in section 6, we discuss the precision of the new models and evaluate their influence over the survival rate of two cylinders during a whole atmospheric reentry, using the spacecraft oriented reentry code ARES.

Wall heat flux and pressure coefficient in the shadow region of space debris

The large number of possible geometries for space debris and the complexity of the flow structures developing on afterbodies in the hypersonic continuous regime make it difficult to apprehend the phenomena occurring in the shadow region of space debris. As a first step for this work, we compared the results from the literature on afterbody hypersonic flows and results from CFD computations. The numerical data comes from a database of Navier-Stokes simulations which was used to develop reduced models for the reentry simulation software DEBRISK [START_REF] Omaly | DEBRISK V3: CNES tool evolutions for re-entry risk analysis[END_REF], developed at CNES. These computations were realized with the CFD solver MISTRAL from R.Tech [START_REF] Tech | [END_REF], for laminar flows in chemical non-equilibrium. The database contains several typical geometries of space debris (cylinders, cubes, flat plates...), and several angles of attack and flight points for each of them. The walls are considered totally catalytic and have a fixed temperature T w = 700 K. Four flow topologies were shown to cause significant levels of C p or heat flux on the walls in the shadow region, and are listed below.

Attached flow is observed in the shadow region when the corresponding wall is parallel or nearly parallel to the incoming flow. Zappa and Reinecke observed that the flow can remain attached for angles between the free flow and the wall up to 28° [START_REF] Zappa | An experimental investigation of base heating on typical Mars entry body shapes[END_REF]. Depending on the geometry and the incoming flow conditions, the maximal heat flux resulting from this flow topology can reach 10 % of the stagnation heat flux [START_REF] Wright | Aerothermal analysis of the project Fire II afterbody flow[END_REF][START_REF] Walpot | Base flow investigation of the Apollo AS-202 command module[END_REF]. Moreover, results from the CNES database have shown that for geometries such as cylinders, boxes or flat plates, a peak heat flux is visible on the lee side near the upwind edge. Using the x axis defined in figure 2, this behavior is visible at x 0.3 m in figure 3 Detached flow with fluid reattachment appears mainly at the base of reentry bodies, where the boundary layer detaches at the trailing edge and reattaches in the wake. As can be seen in figure 4, a recirculation forms behind the base, redirecting the air heated by the recompression at the reattachment point R F towards the afterbody. As a result, one or several heat flux peaks appear on the base wall, depending on the Reynolds number [START_REF] Rom | Laminar heat transfer to a two-dimensional backward facing step from the high-enthalpy supersonic flow in the shock tube[END_REF], where the heat flux can reach as much as 24 % of the stagnation point heat flux [START_REF] Wright | Aerothermal analysis of the project Fire II afterbody flow[END_REF]. This behavior was also observed in the results of the MISTRAL computations, on the aft of cylinders, boxes and flat plates. For instance, figure 5 exhibits a peak of total heat flux on the aft of a cylinder, that reaches 11 % of the stagnation point heat flux.

Detached flow with solid reattachment is mainly observed behind backward facing steps, where the boundary layer detaches at the edge of the step and reattaches on the wall downwind (figure 6). The recompression caused by the reattachment results in a heat flux peak near the reattachment point R S that can reach up to 18 % of the stagnation point heat flux [START_REF] Horvath | Blunt body near wake flow field at Mach 10[END_REF][START_REF] Rom | Laminar heat transfer to a two-dimensional backward facing step from the high-enthalpy supersonic flow in the shock tube[END_REF][START_REF] Hollis | Transition effects on heating in the wake of a blunt body[END_REF].

Shock-shock interactions can happen inside hollow geometries, like truncated cones or pipes. Shock-shock interactions have been classified by Edney [START_REF] Edney | Anomalous heat transfer and pressure distributions on blunt bodies at hypersonic speeds in the presence of an impinging shock[END_REF] into six different types, depending on the angle between the two shocks. Results from the MISTRAL Navier-Stokes computations have shown that type II and III shock-shock interactions can happen in truncated cones and cause high pressure and heat flux peaks on the inside of the cone, up to 60 % of the stagnation point pressure and 50 % of the stagnation point heat flux. These peaks can lead to severe thermo-mechanical stresses on the impacted surfaces, and to strong ablation or fragmentation. The topologies of type II and III interactions between two shocks (C 1 ) and (C 2 ) are drawn in figure 7. This figure also shows in which region of the shock (C 2 ) each type of interaction occurs, in the case where (C 2 ) is the curved shock developing around a sphere, and (C 1 ) is an oblique shock coming from below, for instance from another hypersonic object upstream.

From this primary study, we decided to model the wall C p and heat flux distributions in the shadow region of space debris resulting from attached flows and detached flows with fluid reattachment. Indeed, these two topologies appear on many different geometries, can affect large areas in the shadow region and seem easier to model, compared to solid reattachment and shock-shock interactions. We decided to build the models for the lee side and the aft of cylinders as an initial step and to extend them to boxes and flat plates in future works.

Input parameters and domain of variation

We used the database of Navier-Stokes computations provided by CNES to determine the input parameters of the models, and their domain of variation. We chose to describe the cylinder geometries by their length L and their length-diameter ratio L/D, which has the advantage of being dimensionless. The influence of the edge radius of the cylinder was not visible in the CNES database, since all cylinders had the same edge radius. However, the two modeled phenomena (attached flow on the lee side and detached flow on the aft of the cylinders) involve expansion fans around edges, and the influence of edge radius on detached flow topology was highlighted by O'Byrne [START_REF] O'byrne | Hypersonic laminar boundary layers and near-wake flows[END_REF]. Therefore, it was decided to include the edge radius R e of the cylinders as an input parameter of the models. The domains of variation of the three geometric parameters were chosen as follows:

0 < L < 3 m 0 < L/D < 15 1 mm < R e < 10 mm
The angle of attack α was also chosen as an input parameter of the models, as it strongly influences the wall heat flux and C p levels, both on the lee side and on the aft of the cylinders. For instance, figure 8 shows the distributions of normalized total heat flux on the lee side and on the aft of a cylinder, for values of α between 0 and 45°. For angles α ≥ 15°, the heat flux peak on the lee side of the cylinder is not visible anymore, and the heat flux level is mostly below 5 % of the stagnation heat flux. On the aft of the cylinder, a small peak of heat flux is still visible for α = 15°, but the heat flux level is below 2 % of the stagnation heat flux for all angles of attack α ≥ 15°. Thus, for the models of wall heat flux and C p distributions, we decided to limit the domain of variation to angles of attack between 0 and 15°. The cylinder geometries available in the MISTRAL database all have the same diameter D = 1 m and the same edge radius R e = 6 mm, and 3 different lengths (L = 0.1 m, L = 1 m and L = 3 m). Results are only available for α = 0 or α = 15°, which correspond to the limits of the chosen domain of variation. Therefore, additional numerical computations were required in order to correctly explore the parameters domain of variation.

In order to minimize the number of samples required, we assumed that there was no interaction between the incoming flow conditions and the geometric and attitude parameters L, L/D, R e and α. This assumption is based on correlations from the literature, linking the stagnation point heat flux to the incoming flow characteristics [START_REF] Detra | Addendum to 'Heat transfer to satellite vehicles reentering the atmosphere[END_REF][START_REF] Kemp | Heat transfer to satellite vehicles reentering the atmosphere[END_REF][START_REF] Fay | Theory of stagnation point heat transfer in dissociated air[END_REF][START_REF] Sutton | A general stagnation point convective heating equation for arbitrary gas mixture[END_REF][START_REF] Sagnier | Flow characterization in the ONERA F4 high enthalpy wind tunnel[END_REF]. Thus, we modeled the effects of the incoming flow separately from the other parameters, using existing correlations.

The influence of the remaining four parameters, including their possible interactions, were modeled using POD (Proper Orthogonal Decomposition) and interpolation.

In the following section, we present the bases of the POD and interpolation method, and explain how the entry samples of the model were chosen and computed.

Modeling methodology

Proper Orthogonal Decomposition and interpolation

The basic principle of POD consists in projecting the considered problem on a vector subspace of smaller dimension, chosen to minimize the projection error. Using this method, it is possible to greatly reduce the dimension of the problem, by identifying the most important components of the modeled phenomenon. It is used in several fields of study, such as human face characterization [START_REF] Kirby | Application of the Karhunen-Loève procedure for the characterization of human faces[END_REF] or data compression [START_REF] Andrews | Adaptive data compression[END_REF]. In fluid mechanics, it can be used to predict turbulent structures [START_REF] Sirovich | Turbulence and the dynamics of coherent structures Part 1: coherent structures[END_REF] or even the whole 3D flow around a plane [START_REF] Dupuis | Aerodynamic data predictions for transonic flows via a machine-learning-based surrogate model[END_REF]. In particular, POD was used by Xin et. al. [START_REF] Xin | A reduced order aerothermodynamic modeling framework for hypersonic vehicles based on surrogate and POD[END_REF] to model wall C p and heat flux distributions on a wing as a function of incoming flow speed, altitude and angle of attack.

Considering a distribution w(p) ∈ R N with p a point of the domain of variation Ω and N the dimension of the problem (for instance, the number of points of a CFD mesh), POD consists in building a subspace from R N of dimension n which minimizes the projection error

pro j = Ω w(p) -w(p) 2 dp (1)
where w(p) is the projection of w(p) on the subspace. This subspace is generated by the n eigenvectors (Φ i ) i=1..n associated with the n first eigenvalues λ 1 > λ 2 > . . . > λ n of the autocorrelation matrix

K = Ω w(p)w(p) T dp (2)
When the data w(p) is only known for discrete samples, or "snapshots" (p k ) k=1..N s of the domain Ω, the discrete POD method is used [START_REF] Sirovich | Turbulence and the dynamics of coherent structures Part 1: coherent structures[END_REF]. The autocorrelation matrix K is approached by

K = N s i=1 w(p i )w(p i ) T = SS T (3)
where S = [w(p 1 ) . . . w(p N s )]. Since S S T has the same nonzero eigenvalues as S T S, and since the number of snapshots is often much smaller than the dimension of the problem (N s N), it is generally faster to compute the eigenvectors (Ψ i ) i=1..N s from S T S. Eigenvectors (Φ i ) i=1..N s of K then come from the equation

Φ i = 1 √ λ i S Ψ i (4)
The number n of eigenvectors kept to create the projection subspace is chosen thanks to the eigenvalues (λ i ), that represent the energy corresponding to the associated eigenvector. The parameter n is often chosen so that the ratio of the energy kept over the total energy of the system is greater than a given value, in general

n i=1 λ i N s i=1 λ i > 0.99 (5)
According to this criterion, the three first eigenvectors (Φ i ) i=1..3 were kept to create the models for heat flux and C p distributions in the present study. Indeed, using more eigenvectors did not result in noticeable improvement of the models precision.

Once the eigenvectors of the system have been identified, POD can be used in two ways. The first method consists in projecting the governing equations of the distribution w(p) in the basis formed by the eigenvectors. Approximating w(p) by its projection w(p) in the basis formed by the vectors (Φ i ) i=1..n :

w(p) w(p) = n i=1 a i (p) Φ i (6)
one can replace w(p) by w(p) in the equations of the problem, and obtain a linear system for coefficients (a i ) i=1..n , which can be solved numerically. However, this approach requires to modify the numerical solver of the governing equations of the system, which is not always possible. In some cases, one only has access to a "black box", software or experiment that can determine the solution w(p) for a given point p of the domain Ω, generally at a high cost in time or material. In this case, POD can be used to interpolate the solution w(p * ) for a new point p * from the already known snapshots (w(p i )) i=1..N s [START_REF] Wang | A comparative study of POD interpolation and POD projection methods for fast and accurate prediction of heat transfer problems[END_REF]. The distribution w(p * ) is constructed in the basis of eigenvectors (Φ i ) i=1..n by interpolating the new coefficients (a i (p * )) i=1..n from the coefficients of the snapshots (a i (p k )) i=1..n, k=1..N s .

In practice, it often yields better results to subtract the average of all the snapshots to each snapshot, and to apply POD-interpolation on the variation around the average, since the first mode from POD is often close to the average of all snapshots.

In this study, domain Ω is the previously described domain of variation of the 4 parameters L, L/D, R e and α, and w(p) is the wall C p or heat flux distribution on the lee side or the aft of the corresponding cylinder. The snapshots (w(p i )) i=1..N s are sampled via Navier-Stokes laminar computations, described in section 4.3. The coefficients (a i (p * )) i=1..n of the modeled point are interpolated using kriging, or Gaussian process [START_REF] Krige | A statistical approach to some basic mine valuation problems on the Witwatersrand[END_REF][START_REF] Williams | Gaussian processes for machine learning[END_REF], thanks to the open-source toolbox SMT developed at ONERA in cooperation with University of Michigan (MDOLab), NASA Glenn Research Center and ISAE-SUPAERO [START_REF] Bouhlel | A Python surrogate modeling framework with derivatives[END_REF]. The main idea of kriging is to consider the interpolated quantity y(p) as the realization of a stochastic process conditioned by the sampled data. The nonconditioned process, or prior, can be written as:

Y(p) = µ(p) + Z(p) (7)
The deterministic function µ(p) represents the trend of Y(p), and Z(p) is a Gaussian process centered in 0 whose covariance function can be written as:

Cov(p, p ) = σ 2 r(p, p ) ( 8
)
where r is the spatial correlation function, defined as:

r(p, p ) = exp         d j=1 θ j |p j -p j | δ j         (9)
with d the number of input parameters, θ j > 0 and δ j ∈ [START_REF] Anz-Meador | Orbital debris quarterly news[END_REF][START_REF] Del Portillo | A technical comparison of three low earth orbit satellite constellation systems to provide global broadband[END_REF]. The hyperparameters σ 2 and (θ j ) j=1..d are determined by maximizing the likelihood of Y(p) given the sampled values (y(p i )) i=1..N s . The hyperparameters (δ j ) j=1..d can be determined by a maximum likelihood estimation, but can also be fixed by the user in order to simplify the formulation of the problem and to reduce the computational cost. Depending on the structure of the µ(p) function, three types of kriging algorithms can be defined. Simple kriging assumes that µ(p) is a known constant, ordinary kriging assumes that µ(p) is an unknown constant β 0 , and universal kriging assumes that µ(p) is polynomial:

µ(p) = m-1 j=0 β j f j (p) ( 10
)
where m is the number of basis functions ( f j ) j=0..m-1 , with f 0 = 1, and (β j ) j=0..m-1 are the unknown regression coefficients. For ordinary or universal kriging, the hyperparameters β are determined by maximum likelihood estimation. In the present study, the kriging models were found more accurate when using ordinary kriging and fixing δ j = 1, which corresponds to the exponential correlation function.

Once the hyperparameters of Y(p) have been determined, one can build the conditioned process Ỹ(p) ∼ [Y(p)|Y(p i ) = y(p i ), i = 1..N s ], which is also a Gaussian process. The quantity of interest y(p) is then modeled by the mean E( Ỹ(p)) of Ỹ(p). The main interest of this method is that it gives access to the variance of Ỹ(p) in each point, which is an estimate of the precision of the model between the sample points. In our case, three coefficients (a i (p * )) i=1..3 are modeled, corresponding to the three selected eigenvectors (Φ i ) i=1..3 . These coefficients are modeled separately, and each heat flux or C p distribution model is composed of three kriging sub-models.

Kriging requires an appropriate design of experiments (DOE), i.e. the distribution of the sample points in the domain of variation of the parameters. The design of experiments used in this work is detailed in section 4.2. More details on kriging can be found in Forrester's book on reduced order modeling [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF], among others.

Adaptive design of experiments

As stated before, the effect of altitude h and incoming flow speed U ∞ on wall C p and heat flux distributions is addressed separately, under the hypothesis that there is no interaction with the other parameters. For the remaining parameters L, L/D, R e and α, since the possibility of interactions with one another could not be ruled out by the analysis of MISTRAL computations, the design of experiments (DOE) must cover the whole domain of variation, instead of investigating one parameter at a time. In order to minimize the number of samples needed, we chose to use an adaptive design of experiment: an initial Latin Hypercube Sampling [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF] (LHS) DOE with 12 sample points was constructed for the four parameters (cases n°1 to 12 in table 1). For these 12 first points, the diameter D of the cylinder results from the two other parameters L and L/D which were chosen as entry parameters for the models. However, the primary results exhibited a strong influence of the diameter on the wall heat flux and C p distributions (see section 5). The primary DOE being based on parameter L/D, it resulted in too many cases of small diameters (D < 0.5 m), with several cases of the same diameter. This is clearly visible in figure 9, which shows the projections of the sample points from the initial DOE (in blue) in the planes L -α and L -D. The points are well distributed in the L -α plane, but are too condensed in the region of small values of D in the L -D plane. Thus, 4 samples were added by choosing the parameter D instead of L/D as input parameter. In table 1, the grey boxes show which parameter, from L/D or D, results from the others.

A study of the Sobol' indices [START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF][START_REF] Sobol | On sensitivity estimation for nonlinear mathematical models[END_REF] for the first 12 samples was performed, in order to estimate the importance of each variable. The computation of Sobol' indices usually requires many evaluations of the function of interest. In the present case, we had only 12 evaluations at our disposal, and we performed the evaluation of Sobol' indices using the n surrogate models for coefficients a i built with the 12 samples. Therefore, the computed indices suffer from the same error as the surrogate models, and must be taken with caution. However, we observed that the total indices corresponding to the edge radius R e were below 1 × 10 -2 for all the surrogate models. This made it clear that R e had a negligible influence on the heat flux and C p distributions, both on the lee side and on the aft of cylinders. Thus, for the additional sample points, this parameter was set to R e = 4 mm. This does not invalidate the initial DOE, based on Latin Hypercube Sampling. Indeed, this type of DOE is robust when projected on a subspace, meaning that the samples remain well distributed even when one or several parameters are neglected later in the modeling process. The analysis of the Sobol' indices did not make it possible to rule out any of the remaining three parameters, nor to infirm or confirm the existence of interactions between them.

The choice of the 4 additional samples was made using two different criteria. Samples n°13 and 15 were picked in specific regions of the domain of variation, in order to investigate phenomena that had been observed in the primary results and are detailed in section 5. Sample points n°14 and 16 where chosen using a method of adaptive DOE, based on Integrated Variance-Mean Ratio (IVMR) [START_REF] Liem | Surrogate models and mixtures of experts in aerodynamic performance prediction for mission analysis[END_REF]. According to this method, the new sampling point can be chosen by determining the minimum of the IVMR on the domain of variation Ω:

IVMR(p) = Ω σ 2 p (s) |µ p (s)| + ds ( 11
)
σ 2 p (s) and µ p (s) being respectively the variance and the prediction in point s of the "bogus" model created by taking p as the new sample point, and the prediction of the existing model at point p as the sampled value. In our study, there are as many kriging models as eigenvectors kept for the POD-interpolation. Therefore, we proposed a modified definition of IVMR to take into account all the kriging models. To do so, we considered the mean M (p) of the coefficients (a i (p)) i=1..n weighted by the corresponding eigenvalues (λ i ) i=1..n of the correlation matrix K:

M (p) = n i=1 a i (p) λ i (12)
The weights (λ i ) i=1..n represent the relative importance of the coefficients (a i (p)) i=1..n in the modeled phenomenon. Since the coefficients (a i (p)) i=1..n are modeled using a kriging algorithm, they can be seen as realizations of normal laws of means (µ i,p ) i=1..n and variances (σ 2 i,p ) i=1..n , and their linear combination is also a realization of a normal law, of mean µ Σ,p and of variance σ 2 Σ,p :

                   µ Σ,p = n i=1 µ i,p λ i σ 2 Σ,p = n i=1 σ 2 i,p λ 2 i ( 13
)
Therefore, we defined the IVMR as

IVMR(p) = Ω σ Σ,p 2 (s) |µ Σ,p (s)| ds (14)
Using an optimization algorithm, it is then possible to find the new sampled point p new , which verifies

IVMR(p new ) = min {IVMR(p), p ∈ Ω} (15)
In the present study, we used the surrogate-based gradient-free optimization algorithm SEGOMOE developed at ONERA [START_REF] Bartoli | Adaptive modeling strategy for constrained global optimization with application to aerodynamic wing design[END_REF]. realized for flight points FP1 and FP2 for two cylinder geometries, corresponding to cases n°11 and 15 of table 1. These flight points were chosen from the MISTRAL database, and correspond to a typical debris reentry trajectory. The results for the three flight points were extended to other flow conditions using the Vérant-Sagnier correlation [START_REF] Sagnier | Flow characterization in the ONERA F4 high enthalpy wind tunnel[END_REF], as detailed in section 6.

The final DOE for modeling the effects of geometry and angle of attack only consists of 16 points, which is few considering that there are three entry parameters (since R e is neglected). Indeed, the "rule of thumb" for kriging is that the number of points required is about ten times the dimension of the domain [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF], 30 points in this case. However, the samples of wall C p and heat flux distributions are obtained by expensive 3D chemical non-equilibrium numerical computations, described in section 4.3, and it was not possible to realize this many computations in an affordable time. Moreover, the use of an adaptive DOE algorithm for the choice of samples n°14 and 16 is supposed to reduce the number of samples required to reach a satisfactory level of precision for the wall C p and heat flux models.

FP1

Navier-Stokes computations of sample points

The computations of the sample points were realized with the Eulerian unstructured Navier-Stokes solver CHARME from the software CEDRE developed at ONERA [START_REF] Refloch | CEDRE Software[END_REF][START_REF] Onera Cedre | [END_REF]. According to the flight points considered (velocity, altitude), the flow is assumed in chemical non-equilibrium. The flow reactions are modeled through a 5 species Park model (N 2 , O 2 , N, O, NO) with 17 reactions [START_REF] Park | Nonequilibrium hypersonic aerothermodynamics[END_REF]. The wall is considered totally catalytic with a fixed temperature T w = 700 K. The flow is supposed laminar for the present study, even though Lees' criterion [START_REF] Lees | Hypersonic wakes and trails[END_REF] predicts the transition to turbulent flow on the rear of geometries for flight point FP1. The effect of turbulence on the heat flux distribution is not modeled in the present study. Navier-Stokes equations are solved with a finite-volume discretization and an implicit timestep, since the flow is stationary for the considered angles of attacks. All computations are carried out on full 3D meshes, without making any assumption on the symmetry of the flow. Depending on the cases, AUSM+ or hybrid HLL-HLLC scheme is used. With the latter, HLL scheme and first-order spatial discretization are applied to the shock region, while HLLC and second-order spatial discretization are applied to boundary layer and field. This is made possible by the meshes used in this work, in which the flow field is separated into 3 domains. The boundary layer is meshed with structured hexahedra, the field with unstructured tetrahedra, and the shock with extrusion prisms from the field surface mesh (figure 10). The prisms in the shock region follow the shock geometry and have a maximum height of 2 mm for the biggest geometries. Near the wall, a first cell size of 1 × 10 -6 m high, with an expansion ratio of 1.2, and a maximal length of 15 mm was chosen to properly capture gradients near the wall and thus allow a correct computation of the convective-diffusive heat flux at the wall. For the wake, an initial size of 2.5 mm was chosen based on the work of Barnhardt et. al. [START_REF] Barnhardt | CFD analysis of CUBRC base flow experiments[END_REF], but a mesh convergence study showed that for the considered geometries, a maximum cell size of 10 mm in radial and longitudinal direction yielded similar results, with a maximum difference of 5 % between the heat flux distributions computed with a "coarse" mesh (maximum cell size of 10 mm in the wake) and a "fine" mesh (maximum cell size of 2.5 mm in the wake). Thus, the computations were realized with a maximum cell size in the wake between 2.5 mm and 10 mm depending on the case, in order to keep the total cell number under 50 million. With these cell sizes, a computation required between 40 000 and 250 000 CPU-hours, depending of the cylinder geometry. In the DOE presented in section 4.2, all parameters are varied at the same time. This permits to detect the interactions between the parameters and to reduce the number of required computations, but makes it difficult to give physical explanations to the observed phenomena. Nonetheless, the results of Navier-Stokes computations have highlighted the strong influence of the diameter D of the cylinder on the heat flux distributions on the lee side. This influence can be linked to the Reynolds number based on the cylinder diameter Re D,∞ = U ∞ D/ν ∞ , with U ∞ and ν ∞ the speed and kinematic viscosity of the incoming flow. For cylinders of diameter D ≥ 0.4 m (Re D,∞ ≥ 1.4 × 10 4 ), the heat flux distribution exhibits a peak near the leading edge of the cylinder, as seen in the MISTRAL computations, whereas for cylinders of diameter D ≤ 0.2 m (Re D,∞ ≤ 6.9 × 10 3 ), no peak is visible and the heat flux decreases from the leading edge to the trailing edge (figure 11). The value of D that separates the two behaviors could not be determined, since no case is available with a diameter 0.2 m < D < 0.4 m. This phenomenon seems to be independent of L and α and to be related to the Reynolds number Re D,∞ , which indicates the presence of strong viscous phenomena in the leading edge area. This is confirmed by the results for the other flight points: figure 12 shows the heat flux distributions on the lee side of the cylinder n°15, for the three flight points. It is clearly visible that there is no peak of heat flux for this case at flight point FP3, but that the peak appears at flight point FP2, and is even sharper at flight point FP1. For this geometry (D = 0.15 m), the values of Re D,∞ are 5.2 × 10 3 for flight point FP3 and 1.5 × 10 4 for flight point FP2, which is consistent with the observed behaviors. In conclusion, the heat flux distribution on the lee side of cylinders are monotonous for Re D,∞ ≤ 6.9 × 10 3 and show a peak near the leading edge for Re D,∞ ≥ 1.4 × 10 4 . More computations would be required to determine the limiting value between the two behaviors. similar conclusion can be drawn for the effect of D on the C p distribution. Indeed, for all cases but n°3, the C p distribution shows a peak on the lee side. For instance, figure 13 shows the comparison between the C p distributions along the lee side of cylinders n°3 and n°15. This observation could indicate that for Reynolds number Re D,∞ ≥ 3500, a peak of C p appears near the leading edge of the cylinder. However, only case n°3 has a lower Reynolds number (Re D,∞ = 1700), and other computations would be necessary to confirm this behavior and to determine the limiting value of Re D,∞ between the two types of C p distribution. Still, the computations realized for other flight points seem to corroborate that this phenomenon is also linked to Re D,∞ : as can be seen in figure 14 for case n°11, a higher Reynolds number results in a more pronounced drop in C p near the leading edge followed by a higher peak, as for the heat flux distribution. For cases with Re D,∞ ≥ 1.4 × 10 4 , where both the peaks of heat flux and of C p are visible, the two peaks are not located at the same position on the cylinder. Horvath and Hannemann [START_REF] Horvath | Blunt body near wake flow field at Mach 10[END_REF] observed the same phenomenon in the case of a detached flow with fluid reattachment, but did not provide an explanation for this phenomenon. In the present study, the peak of heat flux is always located upstream of the peak of C p .

The strong influence of Re D,∞ on the heat flux and C p distributions indicates the existence of an interaction between the freestream conditions and the cylinder diameter. This goes against the decision made in section 4 to model the effects of the flight point separately from the other parameters. We detail in section 6 the effect of this interaction on the models for heat flux and C p distributions. he shape of the distributions on the lee side also strongly depends on α: for cases where a peak of heat flux or C p is present, an increase in α will reduce the height of the peak, and move it towards the trailing edge. A bigger value of α also results in a longer distance between the peaks of heat flux and C p , for cases where both peaks are present. The length of the cylinder does not modify the levels and shape of the distributions of heat flux and C p , but it still has to be taken as an input parameter for the models, since it changes the range of abscissa x on which the heat flux and C p are defined. Finally, as stated in section 4, we did not notice any influence of the edge radius R e on the heat flux and C p distributions.

Heat flux and C p distributions on the aft of cylinders

The computed distributions of heat flux and C p on the aft of the cylinders also highlighted the strong influence of the cylinder diameter D, in interaction with angle of attack α. When α = 0, like in case n°1, the flow is fully axisymmetric, and we observe a peak of heat flux and C p at the center of the aft, resulting from the recirculation in the near wake. For nonzero values of α, two profiles of heat flux and C p are observed, depending on the diameter D. For diameters D ≤ 0.4 m (Re D,∞ ≤ 1.4 × 10 4 ), a peak of heat flux and C p is visible, and is situated in the plane of symmetry y = 0, as can be seen in figure 15. For a given diameter, the position of this peak on the z-axis varies with L and α, although no clear relationship was apparent in the available results. The length L of the cylinder also influences the levels of heat flux and C p on the aft: a shorter length results in an aft closer to the front shock, leading to higher levels of temperature and pressure in the near wake, and therefore higher levels of heat flux and C p on the wall.

For diameters D ≥ 1 m (Re D,∞ ≤ 3.5 × 10 4 ) and nonzero values of α, two peaks of heat flux and C p appear on both side of the plane of symmetry y = 0, as visible in figure 16. This is caused by the formation of two counter-rotating vortices in the ake. The formation of several vortices in the wake is often observed for high values of the Reynolds number [START_REF] Wright | Aerothermal analysis of the project Fire II afterbody flow[END_REF][START_REF] Hollis | Transition effects on heating in the wake of a blunt body[END_REF]. In the present cases, the complexity of the flow in the recirculation and of the wall distributions makes it difficult to quantify the effect of L and α on the heat flux and C p . he study of the results for other flight points did not make it possible to quantify the influence of the flight point on the wall distributions, but they confirmed the strong impact of the Reynolds number Re D,∞ on the flow topology in the wake. At flight point FP2, for case n°15 (D = 0.15 m, Re D,∞ = 1.5 × 10 4 ), the heat flux and C p distribution also have one peak in the plane of symmetry y = 0, although this peak is situated higher on the z-axis than at flight point FP3. For case n°11 (D = 0.4 m, Re D,∞ = 3.9 × 10 4 ), even though the Reynolds number is slightly above 3.5 × 10 4 , only one large peak of heat flux and C p is visible at the bottom of the aft (see figure 17), which could indicate that the appearance of secondary vortices in the wake is not only driven by Re D,∞ , but also by L and α. iven the strong physical differences between the two types of heat flux and C p distributions on the aft of cylinders, it seems necessary to build different models for each range of Re D,∞ to accurately represent the two behaviors. Given the number of points in the DOE, it was not possible to build models with a satisfactory accuracy on the aft of cylinders.

Therefore, models were only developed for heat flux and C p distributions on the lee side of cylinders. In section 6, we evaluate the precision of these models, and their impact on a whole reentry trajectory.

Models for heat flux and C p distributions on the lee side of cylinders

As stated in section 4.2, we modeled the effects of the parameters L, D and α separately from the effects of the flight point, in order to reduce the number of required numerical computations. The models for inputs L, D and α at flight point FP3 were built using the POD-interpolation method presented in section 4.1. In section 6.1, we evaluate the precision of these models using the Leave-One-Out cross-validation method [START_REF] Molinaro | Prediction error estimation: a comparison of resampling methods[END_REF]. These models were then extended to other flight points using the Vérant-Sagnier correlation [START_REF] Sagnier | Flow characterization in the ONERA F4 high enthalpy wind tunnel[END_REF] for the heat flux, as detailed in section 6.2. The pressure coefficient being a dimensionless variable, it should be possible to dispense with the flight point effects. Finally, in section 6.3, we demonstrate the impact of the new models on the integrated heat flux received by two different cylinders during their complete atmospheric reentry.

Modeling of the effects of L, D and α

We used the CFD CEDRE results of section 4 and the POD-interpolation method to build 2D models for heat flux and C p on the whole lee side of the cylinders, for input variables L, D and α and flight point FP3. Case n°3 was not used to build the model for C p , since this case is the only one with a Reynolds number Re D,∞ < 3500. This results in a singular C p distribution, which could not be modeled properly by a linear combination of the eigenvectors of the other cases. Cases n°7, 8 and 10 were not used to build the model for total heat flux, because of problems in the distribution obtained with CFD. Even though two types of heat flux distributions were observed, depending on the Reynolds number Re D,∞ , it yielded better results to build one model with all the samples as inputs, rather than one model for each range of Re D,∞ . This is probably due to the fact that we do not have enough samples to build two different models (6 for Re D,∞ ≤ 6.9 × 10 3 and 7 for Re D,∞ ≥ 1.4 × 10 4 ). In this case, the advantage of using more samples for one model beats the disadvantage of modeling two different types of distribution with the same eigenvectors. We also increased artificially the quantity of input data for the kriging models by taking into account the results for the windward side of the cylinders, for cases at incidence. The eigenvectors were computed only from the results on the lee side, and the results from the wind side were projected onto them. We then added the resulting coefficients a i,wind to the input data of our kriging models, by assigning them a negative angle of attack α wind = -α lee . By doing so, we were able to increase the quantity of learning data from our kriging models without modifying the eigenvectors, which represent the physical phenomena on the lee side only. This method allowed us to slightly improve the overall precision of our models without additional computational cost.

To evaluate the precision of the models without performing additional numerical computations, we used the Leave-One-Out cross-validation method [START_REF] Molinaro | Prediction error estimation: a comparison of resampling methods[END_REF]. This method consists in building a sub-model for each sample point of the DOE, using all the other points as input data. The sub-model is then evaluated on the left-over point, resulting in an error value for each sub-model. The mean error of all the sub-models gives an estimate of the error of the real model built using all the sample points. For models with few sample points, as is the case in the present work, this method can lead to overestimate the error for the points on the convex envelope of the DOE, since these points are situated outside of the learning domain of their respective sub-model. However, this method remains an efficient way of estimating the error of a model without having to perform additional measures or computations.

Figure 18 shows the integrated error of the sub-models on the whole lee side of the cylinders, for each case. This error is computed as follows:

= w CFD -w model 2 w CFD 2 ( 16
)
where w CFD is the sampled distribution computed with CEDRE, and w model is the PODinterpolation prediction. The mean error is 22 % for heat flux and 25 % for C p , with a maximal error around 50 %. This comes from the relatively small number of input points, but also from the Leave-One-Out validation method, which results in high error values for points on the convex envelope of the DOE. Still, the precision of most of the sub-models is below 20 %, which tends to indicate a satisfactory precision for the real model. A comparison of the sampled and modeled distributions of heat flux and C p for case n°6 can be seen in figure 19. For this case, the modeling error is 17 % for the heat flux and 16 % for C p , and one can see a good qualitative agreement between the modeled and sampled distributions. As can be seen in figure 20, the levels of C p remain close on most of the lee side for the three flight points. Most of the difference between the flight points concerns the leading edge of the cylinder: as stated in section 5.1, a larger value of Re D,∞ results in more pronounced drop and peak of C p near the leading edge. For case n°15, the position of the peak is also moved slightly downwind, although this phenomenon is not visible for case n°11, and could not be explained. To precisely model the effect of Re D,∞ on the C p distribution near the leading edge, we would need to capture the interactions between the geometrical parameter D and the incoming flow parameters. This would require a larger DOE with incoming flow parameter U ∞ /ν ∞ as additional input parameter, which would lead to many additional computations. Since the effect of Re D,∞ merely concerns a small area near the leading edge of the cylinders, and since the C p distributions remains globally independent of the flight point on the rest of the lee side, we decided to consider C p to be independent from the flight point, and to simply extend our model built for flight point FP3 to all hypersonic flow conditions. oncerning the total heat flux, figure 21 shows the distributions of total heat flux normalized by a reference heat flux taken from Sagnier and Vérant's work [START_REF] Sagnier | Flow characterization in the ONERA F4 high enthalpy wind tunnel[END_REF]:

Q re f = ρ 0.5 ∞ U ∞ ∆H 1.069 (17) ∆H = H tot,∞ -H wall 287.0 × 273.15 (18)
H tot,∞ = 1005.0 * T ∞ + 1 2 U 2 ∞ (19)
where θ = 3000 K. This normalized heat flux is nearly independent of the incoming flow conditions on most of the lee side of the cylinders. A difference between the flight points can still be seen near the leading edge, where the influence of Re D,∞ is the strongest. Like for the C p , modeling this influence would require many more CFD simulations. Moreover, the height and position of the peak of normalized heat flux are relatively unaffected by the incoming flow conditions, as well as the distribution on the rest of the lee side. Therefore, we decided to extend our model for total heat flux at Table 3 shows the relative difference between the C p and normalized heat flux distributions computed with CEDRE for flight points FP1 and FP2, and the distributions computed for flight point FP3. This difference is the error committed by considering that C p and Q tot /Q re f are independent from the flight point. The results show a relatively small error for flight point FP2, which is the closest to flight point atmospheric reentry software ARES from ONERA [START_REF] Prévereaud | Numerical and experimental study of the thermal degradation process during the atmospheric re-entry of a TiAl 6 V 4 tank[END_REF]. This spacecraft oriented code is capable of modeling an atmospheric reentry with 6 degrees of freedom, using the heat flux and C p distributions on the windward areas of an object during its reentry. For the sake of clarity, we chose to compare the integrated heat flux over the whole surface and during the trajectory, with and without the implementation of the new models on the lee side. We modeled the reentry of two hollow cylinders, of respective geometrical parameters L = 1 m, D = 1 m and L = 1.5 m, D = 0.1 m. To determine their masses, the first cylinder is supposed to be made of aluminum with a wall thickness of 3 mm, which results in a mass of 38 kg. The second cylinder is made of titanium TA6V, with a wall thickness of 8 mm, which results in a mass of 15.3 kg.

The initial height is set to 78 km, the initial speed to 7600 m • s -1 and the initial slope angle to -0.1°. These conditions were chosen to reproduce reentry computations with an object-oriented code, in which the reentering object is classically supposed to break into simply shaped fragments at an altitude of 78 km. Moreover, these conditions correspond to a Mach number M ∞ = 26.6 and to Knudsen numbers based on the length Kn L = 3.2 × 10 -3 for the aluminum cylinder and Kn L = 2.1 × 10 -3 for the TA6V cylinder. The initial flight point is thus in hypersonic continuum regime, for which the models have been developed. The integrated heat flux comparison is stopped once M ∞ reaches 5, after which the incoming flow cannot be considered as hypersonic anymore. The initial angle of attack is chosen as α = 0°, which is the incidence that maximizes the heat flux levels on the walls in the shadow area, as illustrated in figure 8. For this incidence, no aerodynamic moment exerts on the cylinder, and the angle of attack remains null during the whole reentry. In an uncontrolled atmospheric reentry, the body is expected to tumble instead of keeping a constant angle of attack. Additional simulations would be required to estimate the contribution of the shadow area to the integrated heat flux in the case of a tumbling piece of debris, but such computations were not carried out in this study. The present example, with α = 0°, represents the extreme case where most of the body walls remain in the shadow area during the reentry. In this case, since the new C p model only applies on the lee side of the cylinders, it has no impact on the aerodynamic forces, and therefore on the trajectories of the cylinders, which can be seen in figure 22. For the heat flux computation, the walls are considered totally catalytic and have a fixed temperature T w = 700 K, since the new model was designed for these boundary conditions. These computations of whole atmospheric reentries with ARES require approximately 0.2 hour of CPU-time and can be carried out on a desktop computer. In comparison, the CFD computation of a single flight point with CEDRE requires around 100 000 hours of CPU-time, including the preliminary computations required to capture the shock's shape. This shows the importance of surrogate models and fast-response codes for the prevision of the on-ground risk posed by orbital debris.

Figure 23 shows the heat flux distribution with and without the new heat flux model for the lee side at the initial point of the reentry. Figure 24 shows the evolution of the integrated heat flux over the whole geometry and all the flight points of the reentry trajectory for the two cylinders, with and without the new heat flux model. A strong difference is visible between the computed heat flux with and without the new model: at the end of the reentry, the integrated heat flux computed for the aluminum cylinder with the new model is 1.5 times higher than without the new models, and 2.5 times higher for the TA6V cylinder. These results show the strong impact that the new models can have on the integrated heat flux, and therefore on the estimated survival rate of debris, especially for long debris with a large shadow area.

Conclusion

To accurately predict the trajectory and the survival rate of reentering space debris, it is necessary to build accurate and computationally affordable models for the wall heat flux and C p distributions. In the present study, we propose such models for the shadow area of the debris, for hypersonic continuous incoming flow. First, we identified four phenomena responsible for relatively high levels of heat flux or C p in the shadow area space debris: attached flow, detached flow with fluid reattachment, detached flow with solid reattachment, and shock-shock interactions. We decided to model the heat flux and C p distributions on the lee side and the aft of cylinders in Figure 24: Total heat flux integrated over the whole surface of two cylinders during the hypersonic continuous phase of their atmospheric reentry, with and without the new heat flux model on the lee side of the cylinders; the solid lines correspond to the total heat flux integrated over time t, and the dashed lines correspond to the total heat flux integrated over altitude h incidence. To do so, we performed 16 CFD simulations to compute the distributions for cylinders with different parameters L, D, R e and α. The parameters for each sample point were chosen using an adaptive DOE. The results of the computations highlighted the influence of the Reynolds number based on the cylinder diameter Re D,∞ on the shape of the heat flux and pressure distributions, both on the lee side and on the aft of the cylinders. Using these results as input data, we built reduced order models for heat flux and C p distributions, by a POD and interpolation method. The models for the lee side performed satisfactorily, and were extended by taking into account the effect of incoming flow on the distributions. Finally, we demonstrated the impact of these models on the integrated heat flux received by two cylinders during their atmospheric reentry. This work confirmed the importance of accurate models in the shadow area to predict the trajectory and survival rate of space debris during their atmospheric reentry. It also highlighted the interest of the POD-interpolation method for building accurate reduced order models while minimizing the number of required CFD simulations.

In future works, additional input C p and heat flux distributions should be computed, to improve the precision of the models on the lee side and to create the models on the aft. The new DOE could include the flight point as an additional parameter, in order to better capture the strong effect of Re D,∞ on the distributions near the leading edge for all freestream conditions. In the long run, more complex phenomena could be modeled, such as shock-shock interactions which can have a strong influence on the survival rates of certain debris. The models could also be improved by adding the effect of turbulence, which can develop in the wake but also on the rear part of long bodies at low altitude [START_REF] Lees | Hypersonic wakes and trails[END_REF]. Turbulence can strongly modify the flow topology and the wall heat flux in the affected areas [6] [14] [41] [42] [START_REF] Holloway | An investigation of heat transfer within regions of separated flow at a mach number of 6[END_REF], but modeling this effect requires to accurately predict the position of the laminar-turbulent transition first, which can be very complex, especially for hypersonic flows.

Acknowledgement

This work was supported by the ONERA and CNES -"Centre National d'Etudes Spatiales" (French Space Agency) through the co-funding of the PhD of V. Drouet. Authors thank CNES for sharing with them their CFD database realized with the MISTRAL CFD code from RTech. This work was granted access to the HPC resources of GENCI under the allocation A0042B10462 -Project mpe0462.

  (a), which shows the distribution of the normalized total heat flux along the lee side of a cylinder of length L = 1 m and diameter D = 1 m. In this case, the peak heat flux reaches 20 % of the stagnation point heat flux. On the rest of the shadow region, the heat flux decreases but remains over 10 % of the stagnation heat flux. The C p distribution, visible in figure 3 (b), exhibits the same behavior.

Figure 2 :

 2 Figure 2: Schematics of a cylinder with notations (the axes x and z lie in the plane of symmetry y = 0 of the cylinder; R e is the edge radius)

Figure 3 :

 3 Figure 3: Total heat flux normalized by the stagnation total heat flux (a) and C p (b) obtained with MISTRAL (CNES database) along the lee side of a cylinder (L = 1 m, D = 1 m) for the freestream conditions α = 0°, h = 70 km, M ∞ = 20, Re ∞ = 3.64 × 10 4 m -1

Figure 4 :

 4 Figure 4: Topology of a detached flow with fluid reattachment behind a cylinder (L = 1 m, D = 1 m) for the freestream conditions α = 0°, h = 70 km, M ∞ = 20, Re ∞ = 3.5 × 10 4 m -1 ; R F is the location of the fluid reattachment

Figure 5 :

 5 Figure 5: Total heat flux normalized by the stagnation total heat flux obtained with MISTRAL (CNES database) along the aft of a cylinder (L = 1 m, D = 1 m) for the freestream conditions α = 0°, h = 70 km, M ∞ = 20, Re ∞ = 3.64 × 10 4 m -1

Figure 6 :

 6 Figure 6: Topology of a detached flow with solid reattachment (from O'Byrne [15]); R S is the location of the solid reattachment

Figure 7 :

 7 Figure 7: Types II and III of shock-shock interactions (C: shock; Σ: shear layer)

  Finally, it was decided to model the influence of the incoming flow for Mach numbers 9 ≤ M ∞ ≤ 20 and for altitudes 40 km ≤ h ≤ 70 km. (a) Lee side (b) Aft

Figure 8 :

 8 Figure 8: Total heat flux normalized by the stagnation total heat flux obtained with MISTRAL (CNES database) on the lee side (a) and on the aft (b) of a cylinder (L = 1 m, D = 1 m) for several values of α, for the freestream conditions h = 70 km, M ∞ = 20, Re ∞ = 3.46 × 10 4 m -1

Figure 9 :

 9 Figure 9: Projections of the 12 samples of the initial DOE (blue) and of the 4 additional samples (red) in planes L -α (left) and L -D (right)

Figure 10 :

 10 Figure 10: Example of hybrid mesh; in orange: shock region, prism elements; in blue: field region, tetrahedral elements; in green: boundary layer and wake, structured hexahedral elements

Figure 11 :

 11 Figure 11: Total heat flux distributions obtained with CEDRE on the lee side of cylinders n°4 (L = 1.5 m, D = 0.2 m, α = 1.5°) and n°13 (L = 3 m, D = 0.4 m, α = 1.5°) for flight point FP3

Figure 12 :

 12 Figure 12: Total heat flux distributions obtained with CEDRE on the lee side of cylinder n°15 (L = 0.6 m, D = 0.15 m) for flight points FP1, FP2 and FP3 and α = 3°A

Figure 13 :

 13 Figure 13: C p distributions obtained with CEDRE on the lee side of cylinders n°3 (L = 0.5 m, D = 0.05 m, α = 1.5°) and n°15 (L = 0.6 m, D = 0.15 m, α = 3°) for flight point FP3 Figure 14: C p distributions obtained with CEDRE on the lee side of cylinder n°11 (L = 2.1 m, D = 0.4 m) for flight points FP1, FP2 and FP3 and α = 5°T

  (a) Total heat flux [kW • m -2 ] (b) C p

Figure 15 :

 15 Figure 15: Total heat flux (a) and C p (b) distributions obtained with CEDRE on the aft of cylinder n°4 (L = 1.5 m, D = 0.2 m) for flight point FP3 and α = 1.5°w

  (a) Total heat flux [kW • m -2 ] (b) C p

Figure 16 :

 16 Figure 16: Total heat flux (a) and C p (b) distributions obtained with CEDRE on the aft of cylinder n°14 (L = 1 m, D = 1.3 m) for flight point FP3 and α = 4°T

Figure 17 :

 17 Figure 17: Total heat flux (a) and C p (b) distributions obtained with CEDRE on the aft of cylinder n°11 (L = 2.1 m, D = 0.4 m) for flight point FP2 and α = 5°G

Figure 18 :

 18 Figure 18: Integrated error of the models for total heat flux (a) and C p (b) on the lee side of cylinders for flight point FP3 (h = 70 km, M ∞ = 20, Re ∞ = 3.46 × 10 4 m -1 ); the red dashed bar symbolizes the mean value

Figure 19 :

 19 Figure 19: Sampled and modeled distributions of total heat flux and C p on the lee side of cylinder n°6 (L = 1.5 m, D = 0.1 m) for flight point FP3 (h = 70 km, M ∞ = 20, Re ∞ = 3.46 × 10 4 m -1 ) and α = 10°C

Hwall 2 ×

 2 wall = c p,t,wall T wall[START_REF] Sutton | A general stagnation point convective heating equation for arbitrary gas mixture[END_REF] c p,t,wall = 1005.0 + 296.0 θ T e θ/T wall e θ/T wall -1 2

Figure 20 :

 20 Figure 20: C p distributions obtained with CEDRE on the lee side of cylinders n°11 (a) and n°15 (b) for flight points FP1, FP2 and FP3

Figure 21 :

 21 Figure 21: Total heat flux distributions obtained with CEDRE and normalized by a reference heat flux on the lee side of cylinders n°11 (a) and n°15 (b) for flight points FP1, FP2 and FP3

Figure 22 :

 22 Figure 22: Atmospheric reentry trajectories computed with ARES for cylinders in aluminum (L = 1 m, D = 1 m, e = 3 mm) and TA6V (L = 1.5 m, D = 0.1 m, e = 8 mm)

Figure 23 :

 23 Figure 23: Heat flux distributions computed with ARES for a hollow cylinder in aluminum (L = 1 m, D = 1 m, e = 3 mm), at α = 0°, h = 70 km and U ∞ = 7000 m • s -1

  (a) L = 1 m, D = 1 m, e = 3 mm, aluminum (b) L = 1.5 m, D = 0.1 m, e = 8 mm, TA6V

table 2 .

 2 Computations were also

	Case	(m)	L/D	D (m)	R e (mm)	α (°)
	n°L
	1	1	1	1	6	0
	2	1	1	1	6	5
	3	0.5	10	0.05	4	1.5
	4	1.5	7.5	0.2	10	1.5
	5	2.8	14	0.2	7	7
	6	1.5	15	0.1	1	10
	7	0.8	4	0.2	3	13
	8	0.3	1.5	0.2	6	4
	9	1	10	0.1	9	11
	10	2.6	6.5	0.4	4	15
	11	2.1	5.25	0.4	1	5
	12	2.4	1.33	1.8	7	7.5
	13	3	7.5	0.4	4	1.5
	14	1	0.8	1.3	4	4
	15	0.6	4	0.15	4	3
	16	1.8	1.2	1.5	4	14

Table 1 :

 1 Parameters of the computed cases; values indicated in the gray boxes are computed from the other parameters (in white boxes)

Table 2 :

 2 Density ρ ∞ (kg • m -3 ) 3.71 × 10 -3 3.64 × 10 -4 9.49 × 10 -5 Freestream characteristics of the flight points of interest

	FP2	FP3

FP3. The error is larger but still acceptable for flight point FP1, for which the effects of Re D,∞ near the leading edge are more important. The larger error for point FP1 may also be due to the fact that the Mach number M ∞ = 9 for this point, which is relatively low in the hypersonic regime. Therefore, the hypersonic effects in the flow may be less important for flight point FP1 than for points FP2 and FP3.

The results also show a larger error for case n°15 compared to case n°11. This probably comes from the fact that cylinder n°15 has a smaller length (L = 0.6 m) compared to n°11 (L = 2.1 m). Therefore, a relatively larger part of the lee side of the cylinder n°15 is affected by the influence of Re D,∞ , which is not modeled here. Using the Leave-One-Out validation method, we compared the C p and heat flux distributions predicted by the models for cases n°11 and 15 and flight points FP1 and FP2, to the distributions computed with CEDRE for the same cases. The computed error, presented in table 4, regroups the modeling error for flight point FP3, as described in section 6.1, and the error resulting from the extension to other flight points, detailed in table 3. As expected, the error is higher for flight point FP1 than for flight point FP2. It ranges between 18 % and 34 %, which is acceptable but could be improved by modeling more accurately the effect of Re D,∞ on the distributions. Given that this effect results from an interaction between the parameter D and the flight point, modeling it would require to modify the DOE from section 4.2 by adding flight point as an input parameter.

Influence of the new models on the full atmospheric reentries of two cylinders

In this section, we demonstrate the influence of our new models on the heat flux received by two cylinders during their atmospheric reentry. To do so, we compare the evolution during the reentry of the total heat flux integrated over the whole geometry with and without the new models for the distributions on the lee side, using the