# Development of new analytical models for pressure and heat flux distribution on space debris afterbodies

Vincent Drouet, Ysolde Prevereaud, Jean-Marc Moschetta

► To cite this version:

Vincent Drouet, Ysolde Prevereaud, Jean-Marc Moschetta. Development of new analytical models for pressure and heat flux distribution on space debris afterbodies. EUCASS 2019, Jul 2019, MADRID, Spain. 10.13009/EUCASS2019-508. hal-02320363

HAL Id: hal-02320363

https://hal.science/hal-02320363

Submitted on 2 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Development of new analytical models for pressure and heat flux distributions on space debris afterbodies

Vincent Drouet*†, Ysolde Prévereaud *, Jean-Marc Moschetta†

*ONERA

ONERA Toulouse Center, 2 avenue Edouard Belin, 31055 Toulouse, FRANCE

vincent.drouet@onera.fr · ysolde.prevereaud@onera.fr

†ISAE

ISAE-SUPAERO, 4 Avenue Edouard Belin, 31055 Toulouse, FRANCE

jean-marc.moschetta@isae-supaero.fr

## Abstract

In order to better predict the risk induced by orbital debris during their atmospheric re-entry, phenomena in the “shadow region” of the debris, i.e. the region that is not directly impinged upon by the freestream flow, were investigated. Four types of flow were shown to cause significant values of heat flux in the shadow region, for continuous hypersonic hyperenthalpic incoming flow conditions : attached flow, detached flow with fluid or solid reattachment, and shock-shock interactions. Models for pressure and heat flux distributions resulting from attached flows on cylinders leeward side were realized using Proper Orthogonal Decomposition interpolation.

## 1. Introduction

Since 1957 and the orbital performance of the satellite Spoutnik-1, the human activity in space has generated a growing number of space debris. Currently, about 15 000 debris larger than 1 cm are in Earth orbit.¹ During the last forty years, 16 000 tons of space debris ranging from ten microns to several meters long have performed a terrestrial atmospheric re-entry. Between 10 and 40% of the reentering mass is estimated to have reached the Earth,² representing a potential threat to ground safety. 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 French Guyana spaceport.

To quickly estimate the on-ground risk posed by a space debris during atmospheric re-entry, it is required to efficiently predict the pressure and heat flux distributions on the debris surface, which are responsible for its ablation, fragmentation and trajectory. The softwares predicting the on-ground risk posed by space debris are classified into two types, "spacecraft-oriented" or "object-oriented". Spacecraft-oriented codes, like the software ARES developed at ONERA, take into account the whole geometry of the debris. Surrogate models are used to predict wall pressure and heat flux distributions on the whole surface of the debris, from which are inferred aerodynamic forces and moments and local mass ablation at each timestep of the atmospheric re-entry. For the continuous hypersonic part of the re-entry, the models used focus on the windward walls, i.e. the walls directly impacted by the freestream flow, and neglect the "shadow region", formed by all the walls of the object that are masked from the freestream flow (fig. 1). On these regions, the pressure coefficient $C_p$ is often set to 0, which results in a very low heat flux. However, some phenomena described in part 2 can cause relatively high levels of $C_p$ and heat flux in the shadow region of space debris. For instance, some studies¹⁵,²⁵ have shown that the heat flux in the shadow region of a re-entry body could reach 20% of the heat flux at the stagnation point. In such cases, neglecting pressure and heat flux distributions in the shadow region could lead to an overestimation of the risk posed by the debris. Furthermore, the attitude of the debris changes during the atmospheric re-entry, sometimes continuously, in a "tumbling" motion. Thus, some walls of the debris can alternatively be in the windward region and in the shadow region, and the heat flux predicted in the shadow region directly impacts the temperature of these walls, and their possible ablation. Object-oriented softwares are based on the hypothesis that during re-entry, the aerothermodynamic constraints will destroy the links between the pieces of the space debris : at a fixed altitude (generally 78 km), the debris will break up into simple geometric shapes, likes
MODELS FOR PRESSURE AND HEAT FLUX ON SPACE DEBRIS AFTERBODIES

cylinders, boxes or sections of sphere (fig. 2). At each timestep of the re-entry, surrogate models associate to each piece a drag coefficient and an ablated mass, averaged over all the possible attitudes of the piece. The risk prediction software DEBRISK developed by CNES, which will be used as a certification tool in the frame of the LOS, is object-oriented. For this type of softwares, underpredicting the wall pressure and heat flux in the shadow region results in false drag coefficients and in a lower integrated heat flux over the whole debris surface. Therefore, the uncertainty on the debris trajectory is increased and the mass ablation is underpredicted, causing the risk posed by the debris to be overestimated.

<figure><img src="image_2.png" /><figcaption>Figure 1: "Shadow region" (hatched zone) on a re-entry body</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Examples of simple shapes composing a satellite⁵</figcaption></figure>

## 2. Wall pressure and heat flux in the shadow region

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 occuring in the shadow region of space debris. In the present study, four categories of phenomena were identified as causing significant levels of pressure or heat flux on the walls in the shadow region, and are listed hereafter. The phenomena were identified using resources from the litterature as well as a database of Navier-Stokes computations for several geometries and incoming flow conditions characteristic of space debris reentries, realized with the CFD codes CEDRE from ONERA and MISTRAL from R.Tech (provided by CNES), realised for laminar flows in chemical non-equilibrium.

**Attached flow** is observed in the shadow region when the corresponding wall is parallel or nearly parrallel to the incoming flow. The flow can remain attached for angles of attack $\alpha$ up to $28^{\circ}27$ and the corresponding heat flux can reach 7% of the stagnation point heat flux.²⁵ Navier-Stokes computations showed that for full cylinders at small angles of attack (less than 15°), a peak heat flux, reaching 15% of the stagnation point heat flux, is visible on the leeward side near the upwind edge (fig. 3, at normalized curvilinear abscissa $s = 0.11$). Moreover, the integrated heat flux on the walls of the leeward side reaches 33% of the heat flux integrated over the whole debris surface. For angles of attack higher than 15°, however, the heat flux levels on the leeward walls are negligible. Still on cylinders, the wall pressure on the leeward side is nearly the same as the wall pressure on the windward side for small angles of attack. For an angle of attack of 15°, the wall pressure on the leeward side can reach up to 45% of the wall pressure on the windward side. Neglecting the pressure coefficient in this region can lead to an error of 180% on the lift coefficient $C_L$ for some configurations, resulting in a poor modelling of the debris trajectory.

**Detached flow with fluid reattachment** appears mainly at the base of re-entry 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,²⁰ where the heat flux can reach as much as 23% 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 (fig. 5). 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.¹³,¹⁵,²⁰ Solid reattachment was also observed with the Navier-Stokes computations on long geometries at medium angle of attack ($\alpha \approx 15^\circ$), where the flow stays attached on most of the leeward side, but detaches near the trailing edge and reattaches at the edge (fig. 6).

**Shock-shock interactions** can happen inside hollow geometries, like truncated cones or pipes (fig. 7). Shock-shock interactions have been classified by Edney¹⁰ into 6 different types, depending on the angle between the two
MODELS FOR PRESSURE AND HEAT FLUX ON SPACE DEBRIS AFTERBODIES

<figure><img src="image_4.png" /><figcaption>Figure 3: Normalized temperature (a) and heat flux distribution (b) around a cylinder (length $L = 3$ m, diameter $D = 0.4$ m, upstream velocity $U_\infty = 5959.4 \cdot s^{-1}$, upstream temperature $T_\infty = 220.1$ K, upstream pressure $P_\infty = 6.02$ Pa, $\alpha = 1.5°$, laminar, chemical non-equilibrium flow) computed with CEDRE</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Topology of a detached flow with fluid reattachment (from Bertin⁴)</figcaption></figure>

shocks. MISTRAL Navier-Stokes computations have shown that type II and III shock-shock interactions can happen in truncated cones and lead to 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.

In addition to these four identified phenomena, several parameters influence the pressure and heat flux distributions in the shadow region of reentering debris. Turbulence, which can appear on the rear of long geometries for the lowest altitudes, tends to increase wall pressure and heat flux levels, especially in detached flows.⁷, ⁸, ¹³⁻¹⁵, ¹⁷, ²¹ Freestream Reynolds number strongly impacts the topology of detached flows, with secondary vortices appearing in the recirculation zones at higher Reynolds numbers.¹³, ²⁵ Reynolds number also impacts wall pressure and heat flux distribution, but its effects are strongly linked to the laminar or turbulent nature of the flow.⁶, ⁸, ¹², ¹³, ¹⁵, ²⁰, ²¹, ²⁵, ²⁷ Geometric parameters such as nose radius $R_n$ and trailing edge radius $R_{te}$ can also influence the wall pressure distribution in detached flows with fluid reattachment, which tends to increase for low values of $R_n$,⁶ and high values of $R_{te}$.¹⁹ Finally, angle of attack $\alpha$ tends to decrease the heat flux levels on walls in areas with detached flow.¹⁵, ²⁷
MODELS FOR PRESSURE AND HEAT FLUX ON SPACE DEBRIS AFTERBODIES

<figure><img src="image_6.png" /><figcaption>Figure 5: Topology of a detached flow with solid reattachment (from O’Byrne¹⁹)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Solid reattachment on the trailing edge of a 2D flat plate ($U_\infty$ = 5959.4 m·s⁻¹, $T_\infty$ = 220.1 K, $P_\infty$ = 6.02 Pa, $\alpha$ = 15°), computed with CEDRE</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Examples of pipe (left) and truncated cone (right) geometries from the MISTRAL database</figcaption></figure>

### 3. Modelling the wall pressure and heat flux distributions resulting from attached flow in the shadow region

#### 3.1 Parameters and domain of variation

The present paper focuses on the modelling of the pressure and the heat flux distributions for walls in the shadow region where the flow remains attached. As stated before, this phenomenon can cause significant levels of pressure and heat flux on long walls nearly parallel to the flow. The analyse of boxes, flat plates and cylinders of the MISTRAL computations dataset permitted to select the parameters of the geometry and the incoming flow that most influence the wall pressure and heat flux distribution. Within this work, model for cylinders have been realised and will be extended to boxes and flat plates in a future work. For cylindrical geometries, the selected parameters are the length *L*, the length-diameter ratio *L/D*, the edge radius *R*<sub>e</sub>, and the angle of attack *α*, with *α* = 0 when the freestream velocity is aligned with the *x*-axis (fig. 3 (a)). The influence of the flow speed *U*<sub>∞</sub> and of the altitude *z* is taken into account by considering several flight points from a typical trajectory of space debris atmospheric re-entry, which are described in subsection 3.3. Depending on the flow speed, the altitude and the length of the debris, it is likely that turbulent flow could develop on the rear of some geometries, but this phenomenon is not addressed in the present paper and will be investigated in the near future. The range of interest for *α* was chosen from 0 to 15°. These values were determined
<figure><img src="image_9.png" /><figcaption>Figure 8: Types II and III of shock-shock interactions (C : shock ; Σ : shear layer)</figcaption></figure>

from the MISTRAL computations, so that the resulting heat flux in the shadow region reaches at least 10% of the stagnation point heat flux. The range for length and diameters where chosen to adequately represent the most frequent geometries of satellite debris. Thus, the length L is chosen between 0.5 m and 3 m, and the diameter is fixed such that the ratio L/D ranges from 1 to 15. Finally, the edge radius R_e ranges from 1 mm to 10 mm.

## 3.2 Modelling method : Proper Orthogonal Decomposition Interpolation

The dataset of MISTRAL computations was realized for a vast number of geometries, attitudes and flight points, and does not include enough cases in the range of the selected parameters, meaning that dedicated numerical computations have to be performed. Considering the number of selected parameters, it was decided to treat the influence of the flight point separately from the other parameters, using correlations based on the literature.²² This leaves 4 parameters to investigate (L, L/D, R_e and α), including their possible interactions with each other. In order to minimize the number of required computations and to correctly sample the whole domain of variation, the surrogate model for wall pressure and heat flux distributions is realized via POD (Proper Orthogonal Decomposition) Interpolation. This method permits to greatly reduce the dimension of a problem. For this reason, it is widely used in aerodynamics, to model wall pressure and heat flux distributions on a wing²⁶ or even full flow structures,⁹,²³ but also in other scientific fields, like image processing.¹⁶

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

$$ \epsilon_{proj} = \int_{\Omega} \|w(p) - \bar{w}(p)\|_2 dp \quad (1) $$

$\bar{w}(p)$ being the projection of $w(p)$ on the subspace. This subspace is generated by the n eigenvectors $(\Phi_i)_{i=1..n}$ associated with the n first eigenvalues $\lambda_1 > \lambda_2 > \dots > \lambda_n$ of the autocorrelation matrix

$$ \hat{\mathbf{K}} = \int_{\Omega} w(p)w(p)^T dp \quad (2) $$

When the data $w(p)$ is only known for discrete samples, or "snapshots" $(p_k)_{k=1..N_s}$ of the domain $\Omega$, the discrete POD
MODELS FOR PRESSURE AND HEAT FLUX ON SPACE DEBRIS AFTERBODIES

method is used.²³ The autocorrelation matrix $\hat{\mathbf{K}}$ is approached by

$$ \mathbf{K} = \sum_{i=1}^{N_s} \mathbf{w}(p_i)\mathbf{w}(p_i)^T = \mathbf{S}\mathbf{S}^T \quad (3) $$

where $S = [\mathbf{w}(p_1) \dots \mathbf{w}(p_{N_s})]$. Since $\mathbf{S}\mathbf{S}^T$ has the same nonzero eigenvalues as $\mathbf{S}^T\mathbf{S}$, and since the number of snapshots is often much smaller than the dimension of the problem ($N_s \ll N$), it is generally faster to compute the eigenvectors ($\mathbf{\Psi}_i$)$_{i=1..N_s}$ from $\mathbf{S}^T\mathbf{S}$. Eigenvectors ($\mathbf{\Phi}_i$)$_{i=1..N_s}$ of $\mathbf{K}$ then come from the equation

$$ \mathbf{\Phi}_i = \frac{1}{\sqrt{\lambda_i}} \mathbf{S} \mathbf{\Psi}_i \quad (4) $$

The number $n$ of eigenvectors kept to create the projection subspace is chosen thanks to the eigenvalues ($\lambda_i$), that represent the energy corresponding to the associated eigenvector. $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

$$ \frac{\sum_{i=1}^{n} \lambda_i}{N_s} > 0.99 \quad (5) \\ \sum_{i=1}^{n} \lambda_i $$

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 $\mathbf{w}(p)$ in the basis formed by the eigenvectors. Approximating $\mathbf{w}(p)$ by its projection $\bar{\mathbf{w}}(p)$ in the basis formed by the vectors ($\mathbf{\Phi}_i$)$_{i=1..n}$:

$$ \mathbf{w}(p) \approx \bar{\mathbf{w}}(p) = \sum_{i=1}^{n} a_i(p) \mathbf{\Phi}_i \quad (6) $$

one can replace $\mathbf{w}(p)$ by $\bar{\mathbf{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 know the governing equations of the system, which is not always possible. In some cases, one only has access to a "black box", a software or experiment that can determine the solution $\mathbf{w}(p)$ for a given point $p$ of the domain $\Omega$, generally at a high cost in time or material. In this case, POD can be used to interpolate the solution $\mathbf{w}(p^*)$ for a new point $p^*$ from the already known snapshots ($\mathbf{w}(p_i)$)$_{i=1..N_s}$. The distribution $\bar{\mathbf{w}}(p^*)$ is constructed in the basis of eigenvectors ($\mathbf{\Phi}_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. During this work, it was noted that applying POD interpolation to the logarithm of the sampled distribution yielded better results in some cases. Indeed, one can then write

$$ \ln(\mathbf{w}(p)) = \sum_{i=1}^{n} a_i(p) \mathbf{\Phi}_i \quad (7) $$

$$ \mathbf{w}(p) = \exp\left(\sum_{i=1}^{n} a_i(p) \mathbf{\Phi}_i\right) = \prod_{i=1}^{n} \exp(a_i(p) \mathbf{\Phi}_i) \quad (8) $$

The modelled distribution $\mathbf{w}(p)$ can now be seen as the deformation of a standard distribution $\exp(a_0(p) \mathbf{\Phi}_0)$, which fitted better the behaviour of the distributions considered in the present paper.

In this study, the domain $\Omega$ is the previously described domain of variation of the 4 parameters $L, L/D, R_e$ and $\alpha$, and $\mathbf{w}(p)$ is the wall pressure or heat flux distribution on the leeward side of the corresponding cylinder. The snapshots ($\mathbf{w}(p_i)$)$_{i=1..N_s}$ are sampled via Navier-Stokes laminar computations, described in subsection 3.4. The coefficients ($a_i(p^*)$)$_{i=1..n}$ of the modelled point are interpolated using kriging, thanks to an opensource toolbox developed at ON-ERA in cooperation with University of Michigan (MDOLab), NASA Glenn Research Center and ISAE-SUPAERO.²⁴ The main idea of kriging is to consider the interpolated quantity as the realization of a stochastic process:

$$ Y(p) = \mu(p) + Z(p) \quad (9) $$
where $Y(p)$ is the random variable associated with the deterministic quantity $y(p)$ being modelled. $\mu(p)$ is the deterministic part of $Y(p)$ and $Z(p)$ is a normal variable with a mean of 0 and with a covariance function depending on the distance between sample points. One of the interests of kriging is that it gives access to the variance of the normal variable $Z(p)$ in each point, which gives an estimate of the precision of the model between the sample points. More details on kriging can be found in Forrester's book on surrogate modelling,¹¹ amongst others.

### 3.3 Design of computational experiment

As stated before, the effect of altitude $Z$ and incoming flow speed $U_\infty$ on wall pressure and heat flux distributions is addressed separately, under the hypothesis that there is no interaction with the other parameters. The computations were made for 3 different flight points, corresponding to a generic debris atmospheric re-entry, and detailed in table 1. For the remaining parameters $L$, $L/D$, $R_e$ and $\alpha$, since the possibility of interactions with one another could not be rolled out by the analysis of MISTRAL computations, the design of experiment (DOE) must cover the whole domain of variation, instead of investigating one parameter at a time. For this reason, and in order to minimize the number of samples needed, a Latin Hypercube Sampling (LHS) DOE was constructed for the four parameters. An initial DOE was realized with 12 points, and 3 points were added afterwards to improve the model. It is to be noted that the sample 14 was chosen using an adaptative DOE algorithm based on Integrated Variance-Mean Ratio (IVMR).¹⁸ According to this method, the new sampling point can be chosen by determining the minimum of the IVMR on the domain of variation $\Omega$:

$$IVMR(p) = \int_{\Omega} \frac{\sigma_p^2(s)}{|\mu_p(s)| + \epsilon} ds \qquad (10)$$

$\sigma_p^2(s)$ and $\mu_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. 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 \in \Omega\} \qquad (11)$$

The resulting DOE is presented in table 2. The computations of the defined cases were realized for flight point A3 only (tab. 1). Computations for flight points A1 and A2 are limited to a few cases. The DOE only consists of 15 points, which is very few considering that there are four entry parameters. Indeed, the "rule of thumb" for kriging is that the number of points required is about ten times the dimension of the domain, 40 points in this case. However, it was not possible to realize this many 3D chemical non-equilibrium numerical computations in an affordable time, and the use of an adaptative DOE algorithm (for choosing point 14) is supposed to reduce the number of points needed. Even with the small number of points available, the models presented in subsection 3.6 perform rather well on the domain of variation of the parameters.

<table><caption>Table 1: Selected flight points</caption><thead><tr><th></th><th>A1</th><th>A2</th><th>A3</th></tr></thead><tbody><tr><td>Altitude Z (km)</td><td>40</td><td>58</td><td>70</td></tr><tr><td>Vitesse U<sub>∞</sub> (m·s<sup>-1</sup>)</td><td>2888.7</td><td>4769.7</td><td>5959.4</td></tr><tr><td>Température T<sub>∞</sub> (K)</td><td>256.26</td><td>250.61</td><td>220.1</td></tr><tr><td>Pression P<sub>∞</sub> (Pa)</td><td>272.72</td><td>26.3</td><td>6.02</td></tr><tr><td>Masse volumique ρ<sub>∞</sub> (kg·m<sup>-3</sup>)</td><td>3.71 × 10<sup>-3</sup></td><td>3.64 × 10<sup>-4</sup></td><td>9.49 × 10<sup>-5</sup></td></tr><tr><td>M<sub>∞</sub></td><td>9</td><td>15</td><td>20</td></tr><tr><td>Re/m (m<sup>-1</sup>)</td><td>6.56 × 10<sup>5</sup></td><td>9.86 × 10<sup>4</sup></td><td>3.46 × 10<sup>4</sup></td></tr><tr><td>Température de paroi T<sub>w</sub> (K)</td><td>700</td><td>700</td><td>700</td></tr><tr><td>Flux de référence Q<sub>ref</sub> (W·m<sup>-2</sup>)</td><td>10 748</td><td>17 732</td><td>18 501</td></tr></tbody></table>

### 3.4 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. According to the flight points considered (velocity, altitude), the
MODELS FOR PRESSURE AND HEAT FLUX ON SPACE DEBRIS AFTERBODIES

<table><caption>Table 2: Parameters of the computed cases</caption><thead><tr><th>Case n°</th><th>L (m)</th><th>L/D</th><th>D (m)</th><th>R<sub>e</sub> (mm)</th><th>α (°)</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>1</td><td>1</td><td>6</td><td>0</td></tr><tr><td>2</td><td>1</td><td>1</td><td>1</td><td>6</td><td>5</td></tr><tr><td>3</td><td>0.5</td><td>10</td><td>0.05</td><td>4</td><td>1.5</td></tr><tr><td>4</td><td>1.5</td><td>7.5</td><td>0.2</td><td>10</td><td>1.5</td></tr><tr><td>5</td><td>2.8</td><td>14</td><td>0.2</td><td>7</td><td>7</td></tr><tr><td>6</td><td>1.5</td><td>15</td><td>0.1</td><td>1</td><td>10</td></tr><tr><td>7</td><td>0.8</td><td>4</td><td>0.2</td><td>3</td><td>13</td></tr><tr><td>8</td><td>0.3</td><td>1.5</td><td>0.2</td><td>6</td><td>4</td></tr><tr><td>9</td><td>1</td><td>10</td><td>0.1</td><td>9</td><td>11</td></tr><tr><td>10</td><td>2.6</td><td>6.5</td><td>0.4</td><td>4</td><td>15</td></tr><tr><td>11</td><td>2.1</td><td>5.25</td><td>0.4</td><td>1</td><td>5</td></tr><tr><td>12</td><td>2.4</td><td>1.33</td><td>1.8</td><td>7</td><td>7.5</td></tr><tr><td>13</td><td>3</td><td>7.5</td><td>0.4</td><td>4</td><td>1.5</td></tr><tr><td>14</td><td>1</td><td>0.8</td><td>1.3</td><td>4</td><td>4</td></tr><tr><td>15</td><td>0.6</td><td>4</td><td>0.15</td><td>4</td><td>3</td></tr></tbody></table>

flow is assumed in chemical non-equilibrium. The flow reactions are modelled through a 5 species Park model (N₂, O₂, N, O, NO) with 17 reactions. 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 it is likely that flow transition should occur on the rear of the longest geometries at the lowest altitudes. This phenomenon will be considered in a coming study, and the current model will be amended to take the effect of turbulence into account.

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. Depending on the cases, AUSM+ or hybrid HLL-HLLC scheme is used. With the latter, HLL scheme and first-order spatial discretisation are applied to the shock region, while HLLC and second-order spatial discretisation 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 hexaedra, the field with unstructured tetraedra, and the shock with extrusion prisms from the field surface mesh (fig. 13). The prisms in the shock region follow the shock geometry and have a maximum height of 2 mm for the biggest geometries. Concerning the boundary layer cells, the first cell is $1 \times 10^{-6}$ m high, with an expansion ratio of 1.2, and a maximal length of 15 mm was chosen. For the wake, an initial size of 2.5 mm was chosen based on literature,³ but the mesh convergence study realized at the beginning of the present work showed that for the considered geometries, a maximum cell size of 10 mm in radial and longitudinal direction yielded similar results. Thus, the computations were realised 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 millions. It should be noticed that the refinement of the mesh in the wake area makes it possible to use the computations also for modelling pressure and heat flux distributions caused by detached flow with fluid reattachment.

### 3.5 Analysis of the computations results

For all the computed cases, the flow remains attached on the leeward side of the cylinder. For all cases but one, the pressure distribution has the same behaviour (see fig. 9) : a pressure drop near the leading edge, caused by the expansion fan that forms at the edge, followed by a pressure peak and a slow decrease. The levels of the drop and peak, as well as their position along the leeward side, depend mostly on *D* and *α*. However, for case 3, the pressure distribution shows no peak and keeps decreasing starting from the leading edge (fig. 10). This is related to the fact that case 3 is the case with the smallest diameter (*D* = 0.05 m).

Two different behaviours also appear for heat flux distribution : for cases with diameters *D* ≤ 0.2 m, the heat flux
MODELS FOR PRESSURE AND HEAT FLUX ON SPACE DEBRIS AFTERBODIES

<figure><img src="image_10.png" /><figcaption>Figure 9: Computed pressure distribution along the middle line of the leeward side for case 13 ($L = 3$ m, $D = 0.4$ m, $R_e = 4$ mm, $\alpha = 1.5^\circ$)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: Computed pressure distribution along the middle line of the leeward side for case 3 ($L = 0.5$ m, $D = 0.05$ m, $R_e = 4$ mm, $\alpha = 1.5^\circ$)</figcaption></figure>

shows no peak and keeps decreasing starting from the leading edge (fig. 11), while for cases with diameters $D \ge 0.4$ m, the heat flux distribution shows a drop followed by a peak and a slow decrease (fig. 12). Again, no clear explanation for this behaviour appears from the flow field analysis. Especially, the phenomenon does not seem to be caused by a difference in molecules recombination at the wall, since the field concentration for each species is nearly the same for all the cases. For cases with $D \ge 0.4$ m, the levels and position of the drop and peak in heat flux distribution seem to depend mostly on $D$ and $\alpha$. This observation was confirmed by the analysis of the Sobol indices of the results, presented in subsection 3.6.

<figure><img src="image_12.png" /><figcaption>Figure 11: Computed heat flux distribution along the middle line of the leeward side for case 9 ($L = 1$ m, $D = 0.1$ m, $R_e = 9$ mm, $\alpha = 11^\circ$)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: Computed heat flux distribution along the middle line of the leeward side for case 11 ($L = 2.1$ m, $D = 0.4$ m, $R_e = 1$ mm, $\alpha = 5^\circ$)</figcaption></figure>

## 3.6 Validation of the method and discussion

Two types of models were developed in the present study. The first models predict the pressure and heat flux distributions on the walls of cylinders along the middle longitudinal line of the leeward side (fig. 14). These models make it easier to observe the effect of each parameter, and can be extended to flat plates and boxes. The second models predict the pressure and heat flux distributions on the whole surface of the leeward side of the cylinders. The models are currently valid for flight point A3, and the computations of the other flight points are still in progress. Model validation was realized using the "Leave One Out" (LOO) method, which consists in creating several models leaving out one different sample point at each time, and evaluating each model on the left-out point. This method requires no set of
MODELS FOR PRESSURE AND HEAT FLUX ON SPACE DEBRIS AFTERBODIES

validation samples, and minimizes the number of samples, but it is less conclusive than using a proper set of validation samples, especially for points at the boundaries of the variation domain. The error between sampled and modelled distribution is computed as follows :

$$ \epsilon = \frac{\|\mathbf{w}_{\text{sampled}}(x) - \mathbf{w}_{\text{modelled}}(x)\|_2}{\|\mathbf{w}_{\text{sampled}}(x)\|_2} \qquad (12) $$

<figure><img src="image_14.png" /><figcaption>Figure 13: Example of hybrid mesh (red : shock region, prism elements ; blue : field region, tetra elements ; green : boundary layer + wake, structured hexa elements)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14: Middle of the longitudinal line of the leeward side where the pressure and heat flux distributions are modelled</figcaption></figure>

The use of a LHS as DOE permits to analyse the sensitivity of the modelled distribution to each parameter and their possible interactions using the Sobol indices. These indices view the response function y and the entry parameters $x_i$ as random variables, and are defined as follow:

$$ S_i = \frac{\mathrm{Var}(\mathbb{E}(Y|X_i))}{\mathrm{Var}(Y)} \qquad (13) $$

$$ S_{ij} = \frac{\mathrm{Var}(\mathbb{E}(Y|X_i, X_j))}{\mathrm{Var}(Y)} - S_i - S_j \qquad (14) $$

$$ S_i^T = S_i + \sum_{i \neq j} S_{ij} + \dots = 1 - \frac{\mathrm{Var}(\mathbb{E}(Y|X_{i-}))}{\mathrm{Var}(Y)} \qquad (15) $$

where $\mathbb{E}(Y|X_i)$ is the expected value of Y conditional on $X_i$, i.e. the function associating a value $x_i$ of the $i^{th}$ parameter to the expected value of the random variable Y when $X_i$ is fixed to $x_i$. $S_i$ is the first order Sobol index of the $i^{th}$ parameter, which measures the main effect of the variable alone, without interactions. $S_{ij}$ is the second order parameter, which represents the effect of interactions between variables $X_i$ and $X_j$. It is possible to define indices of higher degrees, which measure the interactions between three or more variables. $S_i^T$ is the total Sobol index of the $i^{th}$ variable, and measures the total effect of the $i^{th}$ variable, including the interaction with other variables. Each index ranges from 0 to 1, 0 meaning that the variable (or interaction) has no impact on the target function. Sobol indices can be computed easily, e.g. with the SALib.analyze Python package, but require a great number of samples of the target function (about 1000). In the current study, the target quantities are the coefficients of the POD interpolation. Since only 15 samples are available, a surface response for these coefficients is built using kriging, and this surface response is sampled to compute Sobol indices. Therefore, the values of Sobol indices must be interpreted carefully, but give an indication on the importance of each variable. Thus, Sobol indices show that the edge radius $R_e$ of the cylinder has actually very little impact on the wall pressure and heat-flux distribution on the pressure side. This contradicts what was expected, since the expansion zone developing right after the leading edge was thought to play a big role in both pressure and heat flux distribution. Sobol indices also show that D and $\alpha$ are the most important variables for pressure and heat flux distribution, and that L only has a small impact. It was also noted that the ratio L/D, which was chosen to design the DOE and led to several cases with the same D value, was actually not very relevant, and that it would have been better to directly design the DOE directly based on D. Figure 15 show the comparison between the DOE used in the
present work (in blue), based on the parameter L/D, and a LHS based on the parameter D instead (in red), the range for parameters L, Re and α remaining unchainged. It is clear that the DOE used for this work leaves areas of the domain of variation of D unexplored, and that new computations should be added using a LHS based on D.

<figure><img src="image_16.png" /><figcaption>Figure 15: Comparison between the DOE used in the present work, based on L/D (blue) and a 4-D LHS based on D (red), in the plane D-α</figcaption></figure>

### 3.6.1 Models for pressure and heat flux distributions along the middle line of the leeward side of cylinders

The model for pressure distribution along the middle line was developed partly with the POD interpolation method, and partly from an analytic expression. The distribution before the peak is modelled via the POD interpolation method, and the position of the peak along the line is modelled via kriging. The pressure distribution near the trailing edge is modelled by the analytic expression

$$ \frac{P(x) - P_{\infty}}{P_{\infty}} = 9x^{-0.8} \times [2.3(D - 0.4) + 1] \quad (16) $$

As an example, the modelled pressure distribution for case 13 is visible in figure 17, along with the computed distribution for this case. These models have a rather good agreement with the original data, with errors between 1% and 14% except for case 12 where it reaches 37%. This problematic behaviour for case 12 comes from the fact that this case is on the upper limit of the variation domain of parameter D (D = 1.8 m), and is quite far from the case with the second largest diameter (case 14, D = 1.3 m). Using the LOO validation method, this point is outside of the sampled domain, which causes the error to surge and makes it difficult to properly evaluate the precision of the model in this region of the domain of variation. Again, the choice of the L/D ratio to design the LHS revealed problematic, with several cases having the same value for D, and other cases being very far from one another.

Pressure distribution for case 3, however, has no peak and keeps decreasing along the leeward side. Since it is the only distribution with such a behaviour, it is not possible to use POD interpolation with this case. The analytical expression of equation 16 fits well the pressure distribution near the trailing edge (fig. 10, x = 0.5 m), but it is not valid near the leading edge (x = 0), resulting in a global error of 29%. This is probably due to the influence of α on the pressure near the leading edge, which is not taken into account in the analytical expression. A new set of samples is needed to determine the limit value of D for which the pressure distribution has no peak anymore, and to improve the precision of the model for this particular behaviour.

Heat flux distribution along the middle line was entirely modelled using the POD interpolation method, but two different models had to be created, one for cases where D ≤ 0.2 m (cases 3 to 9 and case 15) and the other one for cases where D ≥ 0.4 m (cases 1, 2 and 10 to 14). For the cases with smaller values of D, the heat flux decreases continually between leading and trailing edges. As an example, the modelled heat flux distribution along the leeward side for case 9 (D = 0.1 m) is shown in figure 11, along with the computed heat flux distribution for this case. This model fits rather well the sampled distributions, with errors between 9% and 28%, all the more so considering that it was designed
MODELS FOR PRESSURE AND HEAT FLUX ON SPACE DEBRIS AFTERBODIES

<figure><img src="image_17.png" /><figcaption>Figure 16: Sampled and modelled pressure distribution along the middle line of the leeward side for case 13 ($L = 3$ m, $D = 0.4$ m, $R_e = 4$ mm, $\alpha = 1.5^\circ$)</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 17: Sampled and modelled pressure distribution along the middle line of the leeward side for case 3 ($L = 0.5$ m, $D = 0.05$ m, $R_e = 4$ mm, $\alpha = 1.5^\circ$)</figcaption></figure>

from only 8 samples, but it still could be improved with new samples. The second model is designed for cylinders with larger diameters ($D \ge 0.4$ m)), and an example can be seen in figure 19 for case 11. This model does not fit as well the sampled data, with errors between 25% and 29%, except for case 12 where it reaches 120%. Again, case 12 is problematic for the LOO validation method, for the same reason as the one explained previously. It is to be noted that the model for larger diameters was designed with only 7 points, 3 of which having the same diameter $D = 0.4$ m. New Navier-Stokes computations are required to improve the model for the larger diameters, and to determine precisely the limit value of $D$ between the two models.

<figure><img src="image_19.png" /><figcaption>Figure 18: Sampled and modelled heat flux distribution along the middle line of the leeward side for case 9 ($L = 1$ m, $D = 0.1$ m, $R_e = 9$ mm, $\alpha = 11^\circ$)</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 19: Sampled and modelled heat flux distribution along the middle line of the leeward side for case 11 ($L = 2.1$ m, $D = 0.4$ m, $R_e = 1$ mm, $\alpha = 5^\circ$)</figcaption></figure>
### 3.6.2 Models for pressure and heat flux distributions on the whole leeward side of cylinders

The POD-interpolation method also permits to model the whole 2D-distributions of pressure and heat flux on the leeward side of cylinders. Figure 20 and 21 show the sampled and modelled distributions for the pressure and heat flux on the leeward side for case 6. On both figures, the horizontal axis shows the distance along the length of the cylinder, and the vertical axis shows the curvilinear distance around its circumference. Only one half of the cylinder is shown, since all cases are symmetrical about thexz plane. The figures show a good agreement between the modelled and sampled distributions. On all the cases, using the LOO validation method, the error ranges from 6% to 35% for the pressure, and from 6% to 27% for the heat flux, except for cases 12 and 14, where the error on heat flux reaches 55% and 45%. The bigger error on these two cases is caused by the lack of sample points with diameters larger than 1 m. These results show the validity of the POD-interpolation method for the modelling of heat flux and pressure distributions in the shadow region of space debris, but new sample points must be added to improve the models.

<figure><img src="image_21.png" /><figcaption>Figure 20: Sampled and modelled pressure distribution on the leeward side for case 6 ($L$ = 1.5 m, $D$ = 0.1 m, $R_e$ = 1 mm, $\alpha$ = 10°)</figcaption></figure>

## Conclusions

In order to develop surrogate models for wall pressure and heat flux distributions in the shadow region of orbital debris, a dataset of CFD computations for several debris geometries, attitudes and hypersonic flow conditions was analysed and compared with results from the literature. Four phenomena were identified to cause significant levels of heat flux and pressure on the walls of the shadow region: attached flows, detached flows with fluid or solid reattachment, and shock-shock interactions. Models for pressure and heat flux distributions caused by attached flows on the middle line of the leeward side of cylinders were developed using POD interpolation, with four entry parameters ($L, L/D, R_e$ and $\alpha$). 15 samples points were obtained with reactive laminar Naviers-Stokes computations, using the solver CHARME from CEDRE. The model for pressure distribution combines the POD interpolation method with an analytic expression, yielding generally good results with errors below 14%, except for extreme values of $D$. Two models for the heat flux distribution were developed, each with a different domain of validity: $D \le 0.2$ m and $D \ge 0.4$ m. The model for smaller diameters fits rather well the sampled data, with errors between 9% and 30%, even though it was designed from only 8 sampled points. The model for larger diameters does not fit as well the sampled data, with errors between 25% and 29%, except for one case where it reaches 120%. Models were also developed for the pressure and heat flux distributions on the whole leeward side of cylinders, with errors ranging from 6% to 35% for the pressure, and from 6% to 55% for the heat flux. It was shown that the POD-interpolation method can be used to model pressure and heat flux distributions in the shadow region of space debris, even though new samples are required to improve the models. These samples can be added by the adaptive DOE method, using the parameter $D$ instead of $L/D$, and removing the parameter $R_e$, which was shown to have a negligible impact on pressure and heat flux distribution. The effect of
MODELS FOR PRESSURE AND HEAT FLUX ON SPACE DEBRIS AFTERBODIES

<figure><img src="image_22.png" /><figcaption>Figure 21: Sampled and modelled heat flux distribution on the leeward side for case 6 ($L$ = 1.5 m, $D$ = 0.1 m, $R_e$ = 1 mm, $\alpha$ = 10°)</figcaption></figure>

turbulence must be also be added to the model for the lower altitudes and the largest debris length. Finally, the models need to be extended to geometries like flat plates and boxes.

## Acknowledgements

The work presented here was supported by a grant from CNES and ONERA, and the CFD simulations were realised thanks to computing resources provided by GENCI. The authors thank Julien Annaloro (CNES) for the provision of the MISTRAL CFD database, and for his help on geometries representative of space debris. The authors thank Nathalie Bartoli and Sylvain Dubreuil (DTIS, ONERA) for their help on POD, kriging and adaptative DOE.

## References

[1] Orbital debris quarterly news, May 2019.

[2] William Ailor, Lyndon Bonaparte, and Alice F. Shelton. A strategy for reducing hazards from reentry debris. In *AIAA Atmospheric Flight Mechanics Conference and Exhibit*, 2006.

[3] Michael Barnhardt, Graham V. Candler, and Matthew MacLean. Cfd analysis of cubrc base flow experiments. In *48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition*, 2010.

[4] John J. Bertin. *Hypersonic Aerothermodynamics*. AIAA, 1994.

[5] Bitlanders. What's a satellite ?, 2014.

[6] Bruce M. Bulmer. Study of base pressure in laminar hypersonic flow : Re-entry flight measurements. *AIAA Journal*, 13(10):1340-1348, 1975.

[7] Dean R. Chapman. A theoretical analysis of heat transfer in regions of separated flow. Technical report, Ames Aeronautical Laboratory, Moffett Field, Calif., 1956.

[8] Dean R. Chapman, Donald M. Kuehn, and Howard K. Larson. Preliminary report on a study of separated flows in supersonic and subsonic streams. Technical report, NACA, 1956.

[9] Romain Dupuis, Jean-Christophe Jouhaud, and Pierre Sagaut. Aerodynamic data predictions for transonic flows via a machine-learning-based surrogate model. In *2018 AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference*, 2018.
MODELS FOR PRESSURE AND HEAT FLUX ON SPACE DEBRIS AFTERBODIES

[10] B. Edney. Anomalous heat transfer and pressure distributions on blunt bodies at hypersonic speeds in the presence of an impinging shock. Technical report, 1968.

[11] Alexander I. J. Forrester, András Sóbester, and Andy J. Kean. *Engineering Design via Surrogate Modelling: A Practical Guide*. 2008.

[12] A. B. Gorshkov. Laminar near wake of a sharp wedge in hypersonic perfect-gas flow. *Fluid Dynamics*, 45(1):126–133, 2010.

[13] Brian R. Hollis and John N. Perkins. Transition effects on heating in the wake of a blunt body. In *32nd AIAA Thermophysics Conference*, 1997.

[14] Paul F. Holloway, James R. Sterett, and Helen S. Creekmore. An investigation of heat transfer within regions of separated flow at a mach number of 6.0. Technical report, NASA, 1965.

[15] Thomas Horvath and Klaus Hannemann. Blunt body near wake flow field at mach 10. Technical report, NASA, 1997.

[16] M. Kirby and L. Sirovich. Application of the karhunen-loève procedure for the characterization of human faces. *IEEE Transactions on Pattern Analysis and Machine Intelligence*, 12(1):103–108, January 1990.

[17] Robert L. Kruse. Transition and flow reattachment behind an apollo-like body at mach numbers to 9. Technical report, NASA, 1968.

[18] Rhea P. Liem and Joaquim R. R. A. Martins. Surrogate models and mixtures of experts in aerodynamic performance prediction for mission analysis. In *15th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference*, 2014.

[19] Sean O'Byrne. *Hypersonic Laminar Boundary Layers and Near-Wake Flows*. PhD thesis, The Australian National University, 2001.

[20] Josef Rom and Arnan Seginer. Laminar heat transfer to a two-dimensional backward facing step from the high-enthalpy supersonic flow in the shock tube. *AIAA Journal*, 2(2):251–255, 1964.

[21] Josef Rom, Arnan Seginer, Rimon Ariely, and Michael Green. Heat transfer in separated regions in supersonic and hypersonic flows. In *The Eight Congress of the International Council of the Aeronautical Science*, 1972.

[22] Philippe Sagnier and Jean-Luc Vérant. Flow characterization in the onera f4 high enthalpy wind tunnel. *AIAA Journal*, 36(4):522–531, April 1998.

[23] Lawrence Sirovich. Turbulence and the dynamics of coherent structures part 1: Coherent strucutres. *Quarterly of Applied Mathematics*, 1987.

[24] (MDOLab) University of Michigan, Nasa Glenn Research Center, ISAE-SUPAERO, and ONERA. Surrogate modeling toolbox, 2017.

[25] Michael Wright, Mark Loomis, and Periklis Papadopoulos. Aerothermal analysis of the project fire ii afterbody flow. *Journal of Thermophysics and Heat Transfer*, 17(2):240–249, 2003.

[26] Chen Xin, Liu Li, Long Teng, and Yue Zhenjiang. A reduced order aerothermodynamic modeling framework for hypersonic vehicles based on surrogate and pod. *Chinese Journal of Aeronautics*, 2015.

[27] O. L. Zappa and W. G. Reinecke. An experimental investigation of base heating on typical mars entry body shapes. In *AIAA 7th Thermophysics Conference*, 1972.