# Pampero HiFi: improvements in aerodynamic modeling of a spacecraft-oriented destructive re-entry code

Valentin Ledermann, Iko Midani, Eddy Constant, Pierre Van Hauwert, Martin Spel, Laurent Stainier, Julien Annaloro, Stéphane Galera

► To cite this version:

Valentin Ledermann, Iko Midani, Eddy Constant, Pierre Van Hauwert, Martin Spel, et al., Pampero HiFi: improvements in aerodynamic modeling of a spacecraft-oriented destructive re-entry code. CEAS Space Journal, 2025, 10.1007/s12567-025-00655-2, hal-05508973

HAL Id: hal-05508973

https://hal.science/hal-05508973v1

Submitted on 13 Feb 2026

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 CC BY 4.0 - Attribution - International License
# Pampero HiFi: improvements in aerodynamic modeling of a spacecraft-oriented destructive re-entry code

Valentin Ledermann<fnref n="1" />,<fnref n="3" /> · Iko Midani<fnref n="1" /> · Eddy Constant<fnref n="1" /> · Pierre Van Hauwert<fnref n="2" /> · Martin Spel<fnref n="1" /> · Laurent Stainier<fnref n="3" /> · Julien Annaloro<fnref n="4" /> · Stéphane Galera<fnref n="4" />

Received: 26 December 2024 / Revised: 29 July 2025 / Accepted: 30 July 2025
© The Author(s) 2025

## Abstract

This work presents the BLIZZARD CFD module developed in the scope of the implementation of the high-fidelity branch of the PAMPERO spacecraft-oriented re-entry tool. First, the Riemann solver, the immersed boundary method, and the adaptive mesh refinement of BLIZZARD are briefly detailed. The mass loss generated by the chosen immersed boundary method is investigated on a sphere and a satellite case. Then, the results of the simulation of static hypersonic wind tunnel test cases are presented. A part of a cube's trajectory computed by the modified Newton theory (method used in spacecraft-oriented tools) is computed using the pressure loads of BLIZZARD. The impact of the fluid modeling on the stability of the object is shown. Finally, the surface pressure distribution and the aerodynamic coefficients obtained with a body-fitted CFD reference code, the modified Newtonian theory and BLIZZARD on a satellite case are compared. BLIZZARD shows major improvement compared to the modified Newtonian theory.

**Keywords** Atmospheric re-entry · Immersed boundary methods · Adaptive mesh refinement · Design for Demise · Spacecraft-oriented destructive re-entry tool · PAMPERO · BLIZZARD

## Acronyms

<table><tr><td>AMR</td><td>Adaptative mesh refinement</td></tr><tr><td>CFD</td><td>Computational fluid dynamics</td></tr><tr><td>CNES</td><td>Centre national d'études spatiales (French Space Agency)</td></tr><tr><td>DOF</td><td>Degree of freedom</td></tr><tr><td>DSMC</td><td>Direct simulations Monte Carlo</td></tr><tr><td>FEM</td><td>Finite element methods</td></tr><tr><td>FSI</td><td>Fluid-structure interactions</td></tr><tr><td>GMRES</td><td>Generalized minimal RESidual.</td></tr><tr><td>GPFM</td><td>Ghost point forcing method</td></tr><tr><td>IBM</td><td>Immersed boundary methods</td></tr><tr><td>MUSCL</td><td>Monotonic upstream-centered scheme for conservation laws</td></tr><tr><td>VKI</td><td>Von Karman Institute for fluid dynamics</td></tr></table>

## 1 Introduction

For many years, the modeling of atmospheric re-entry of space debris using computational fluid dynamics (CFD) has been a major challenge due to its complexity and practical industrial applications. Ablation phenomena are caused by the high thermal loads involved in re-entry, so using a

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

Valentin Ledermann
valentin.ledermann@rtech-engineering.com

Iko Midani
iko.midani@rtech-engineering.com

Eddy Constant
eddy constant@rtech-engineering.com

Pierre Van Hauwert
pierre@rtech-engineering.nl

Martin Spel
martin.spel@rtech-engineering.com

Laurent Stainier
laurent.stainier@ec-nantes.fr

Julien Annaloro
julien.annaloro@cnes.fr

Stéphane Galera
stephane.galera@cnes.fr

<footnote n="1">¹ R.Tech, 10 Impasse Jean Mermoz, Verniolle 09340, France</footnote>

<footnote n="2">² R.Tech Engineering B.V, Baarsjesweg 243,
Amsterdam 1058AA, The Netherlands</footnote>

<footnote n="3">³ Central Nantes, GEM, 1 rue de la Noë, Nantes 44000, France</footnote>

<footnote n="4">⁴ Centre National d'Etudes Spatiales (CNES), 18 Avenue
Edouard Belin, Toulouse 31400, France</footnote>
"high-fidelity" CFD method to accurately measure forces and heat flux is difficult and time-consuming as it requires manual mesh engineering for each point of the trajectory. Consequently, methods based on the modified Newtonian theory are used to model aerodynamic forces, incorporating statistical correlations from CFD and experiments to estimate heat fluxes, sacrificing accuracy for CPU and engineering time [1]. Moreover, these methods, although extremely efficient in computational time, cannot represent the complex phenomena that occur due to the 3D nature of the physics such as shock interactions. To fill the gap between approximate and high-fidelity methods, manufacturers increasingly demand accurate modeling of aero-thermodynamic phenomena throughout a trajectory with realistic computational time and engineering costs for an industrial application. This can be achieved thanks to the improvement of CPU resources in recent years, which has made solving the fluid in an atmospheric re-entry simulation possible. Consequently, fluid-structure coupling strategies have been integrated into new type of re-entry tools, called "Multi-fidelity" tools [2, 3].

The CFD module **BLIZZARD**, developed in the scope of the implementation of the **PAMPERO HiFi** high-fidelity re-entry tool, solves the Euler equations with the perfect gas assumption and uses an automatic octree Cartesian grid generation method to solve the fluid in a re-entry simulation. Automatic grid generation makes CFD simulations on geometries of any kind possible, where the use of body-fitted meshes is hard to set up. The major drawback of this type of approach is the low quality of the near-wall modeling due to the nature of Cartesian meshes. Many methods have emerged to overcome this problem. Among them are the popular overset methods and immersed boundary methods (IBM). **BLIZZARD** uses an IBM, a type of method that was originally used to model fluid flows around curved and/or moving bodies. IBMs were introduced by Peskin and have been since applied to solve fluid flow problems for many applications [4–6].

Several approaches of IBMs have been developed to be used for specific applications. IBM approaches can be classified into two categories: the diffuse-interface methods and the sharp-interface methods. In the diffuse interface IBM, the boundary is spread by distributing forcing and source terms [7, 8] to the surrounding Cartesian grid cells using a window function. Imposing the Dirichlet boundary conditions accurately is simple using this approach, whereas the imposition of the Neumann boundary condition has proven to be challenging. That is why the diffuse interface is usually implemented in incompressible flow solvers. Examples of applications of this method for compressible solvers can be found in the following references: [8–12].

<figure><img src="image_3.png" /><figcaption>**Fig. 1** Automatic mesh generation at initialization on satellite shape</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Mesh obtained after the refinement process</figcaption></figure>

Alternatively, the sharp-interface methods track the boundaries in fluid simulations. Several methods fall under this category, such as the cut-cell method [13, 14], the Cartesian IBM [15], and the ghost-cell method [16, 17]. Among these, the cut-cell method provides the clearest interface by dividing grid cells into sub-cells for the solid and fluid phases. This enables strict enforcement of mass, momentum, and energy conservation. However, the cut-cell method presents challenges in simulating complex body problems due to the intricate cell reshaping process. In the Cartesian IBM, points within the solid body that have at least one neighboring point in the fluid are known as immersed boundaries or forcing points. By reconstructing flow variables at these points and imposing boundary conditions, the accuracy of the flow field can be ensured. The ghost-cell method operates similarly to the Cartesian IBM, but with ghost points representing solid points adjacent to the fluid domain. The reconstruction procedure occurs at these ghost points. While
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Illustration of the generation of the image points. White triangles are fluid cells and black triangles are solid cells, white circles are the image points, green triangles are the surrounding cells used to create the interpolation support for the image points illustrated by green circles, and ghost cells are illustrated by red circles</figcaption></figure>

the Cartesian IBM and the ghost-cell method offer less accuracy ([18, 19]) compared to the cut-cell method due to their implicit representation of the solid boundary, they have simpler point recognition procedures and do not require flux computations around the immersed boundary. Given that the BLIZZARD solver aims to be functional on industrial complex geometries, the ghost-cell method has been selected.

The BLIZZARD fluid solver is loosely coupled with both a six degrees of freedom (6-DOF) solver and a finite element methods (FEM) solver to compute the fluid-structure interactions. While the 6-DOF solver allows a quick estimation of the global rigid body motion, the FEM solver computes the local displacement and deformation of the geometry, as well as the local internal stress distribution. Fluid-structure interactions (FSI) phenomena such as wake effect or shock "surfing" are important for a better evaluation of debris trajectory. The local stress field can then be used to assess the mechanical fragmentation events during re-entry. Details on the coupling scheme and structural governing equations can be found in [20, 21].

The latest capacities of the BLIZZARD fluid solver in the assessment of pressure loads in the framework of industrial atmospheric re-entry are presented in this article. First, the BLIZZARD solver and the sharp interface IBM are briefly detailed. The mass conservation of the implemented ghost-cell method, along with the dependancy of the mass loss on the mesh size, is then studied on a sphere and a satellite. Numerical simulations of the Von Karman Institute for Fluid Dynamics (VKI)'s hypersonic wind tunnel tests carried out with BLIZZARD are compared to a CFD reference code and experimental data to validate the models. The final section presents two application cases: simulations performed by BLIZZARD and the 6-DOF solver of part of a cube's trajectory computed initially by the spacecraft-oriented tool PAMPERO [22] and simulations of the JASON satellite under re-entry conditions, where the BLIZZARD mesh convergence is verified and the pressure loads and aerodynamic coefficients using BLIZZARD, MISTRAL-CFD [23, 24] and the modified Newtonian theory are compared.

## 2 The BLIZZARD solver

In the continuum regime, the pressure loads in re-entry tools are computed using an analytical model (Modified Newtonian Theory), and heat loads are estimated based on statistical correlations. While these methods are computationally efficient, they cannot accurately represent complex physics such as shock interactions even for simple geometry. To improve the estimation of pressure loads the CNES and R.Tech have developed a simplified CFD solver, called BLIZZARD, which solves the Euler equations on the automatically generated Cartesian grid. The results from the CFD solver are also used to obtain a better estimation of the heat loads thanks to a model based on machine learning that is currently under development [25].

The choice of an Euler solver instead of a Navier-Stokes solver is made for two main reasons. The first is the need for a proper boundary layer grid to be created automatically and to be robust, which is a very challenging task. In the regime of interest, the hypersonic continuum, the shear stresses have

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Left: Slice of the MISTRAL mesh on a satellite case. Right: Wall mesh used for IXV re-entry rebuild</figcaption></figure>
only a small contribution to forces and moments acting on the fragments. The second reason is the increased computation time of Navier-Stokes due to the boundary layer clustering, leading to a higher cell count. Thus, the choices of an Euler solver for the pressure distribution and of a machine learning model for the heat fluxes were selected. In the following, numerical methods are briefly exposed, and more details on the BLIZZARD models can be found in [20] and [25].

**Numerical methods** The BLIZZARD module solves the Euler equations considering a perfect gas. For spatial discretization, a finite-volume approach is implemented. For steady computations, a first-order upwind scheme is used, while a third-order Runge-Kutta scheme is employed for unsteady computation. Inviscid flux at each face is estimated using a hybrid VanLeer-AUSM scheme [25].

**Adaptative mesh refinement** BLIZZARD employs a linear octree mesh generation technique using the Bit-pit library [26] to automatically generate a non-uniform Cartesian mesh, enabling efficient execution. This process is based on a patch-based method that identifies fluid and boundary cells during initialization. To detect which cells should be identified as a boundary cells, a surface detection algorithm based on a ray-tracing method powered by the Embree library [27] is used. This algorithm accurately identifies surface nodes that intersect the fluid domain and can handle any type of geometry. Figure 1 shows the resulting initial Cartesian mesh generated on a satellite shape.

The solver uses adaptative mesh refinement (AMR) to improve the accuracy where needed so that more results can be obtained. The adaptation criterion is based on a multi-variable gradient. Consequently, some cells are tagged for refinement. The gradient-based sensor triggers the mesh adaptation mainly in shocks region. Figure 3 illustrates the mesh obtained at the end of the simulation on a satellite shape (Fig. 2).

**Immersed boundary method** One of the primary limitations of this type of approach arises from the lower quality of the near-wall modeling due to the nature of Cartesian meshes. Numerous methods have been developed to address this issue. BLIZZARD uses a ghost point forcing method (GPFM), an IBM that involves the direct imposition of primitive flow variables on the cells surrounding the geometry. Although some formulations impose the flow variables on the cell identified as inside the geometry, the present method is inspired by previous research [28], in which the flow variables are imposed on the solid cells, which become ghost cells.

To impose the boundary condition, an image point is created for each ghost cell by projecting the cell center along the surface normal, see Fig. 3. The flow field variables for each image point are then interpolated from the surrounding cells using an interpolation method based on inverse distance weighting [29]. Once the flow variables are interpolated at the

<figure><img src="image_7.png" /><figcaption>Fig. 5 Surface mesh of the sphere</figcaption></figure>

<table><caption>**Table 1** Free-stream conditions, sphere case</caption><tr><td>M<sub>∞</sub></td><td>6.97</td></tr><tr><td>ρ<sub>∞</sub> [Pa]</td><td>7.835e<sup>-3</sup></td></tr><tr><td>u<sub>∞</sub> [m.s<sup>-1</sup>]</td><td>1046</td></tr><tr><td>T<sub>∞</sub> [K]</td><td>56</td></tr></table>

image point, the boundary condition is imposed by applying the interpolated value to the corresponding ghost cell. In our case, an adiabatic wall condition, a zero-pressure gradient, and a slip condition for the velocity are applied. The flow variables imposed on the ghost cells read:

$$u_g = 2(\dot{d} + u_{ip}^t) - u_{ip}, \quad (1)$$

$$T_g = T_{ip} - \Delta l \left( \frac{\partial T}{\partial n} \right)_b, \quad (2)$$

$$p_g = p_{ip} - \Delta l \left( \frac{\partial p}{\partial n} \right)_b, \quad (3)$$

where subscripts g, ip and b denote variable evaluated at ghost cell, image points and body wall, respectively. $u_{ip}^t = u_{ip} - u_{ip} \cdot n$ is the tangential fluid velocity and $\dot{d}$ is the local solid velocity. For an adiabatic wall $(\partial T / \partial n)_b = 0$ for both steady and moving solids, the pressure gradient is expressed by:
<figure><img src="image_8.png" /><figcaption>**Table 2** Free-stream conditions, instrumented hollow hemisphere</figcaption></figure>

<table><tr><td>M</td><td>11.5</td></tr><tr><td>P[Pa]</td><td>360</td></tr><tr><td>T[K]</td><td>81.4</td></tr><tr><td>u[m/s]</td><td>2124</td></tr><tr><td>AoA[°]</td><td>90.0</td></tr></table>

### 3 Reference code and presentation of the cases

This section presents the numerical methods and the meshing strategy of the reference CFD solver MISTRAL used in this paper and the test cases analyzed in this paper. For each one of them, the geometry model is described, and the inflow/initial conditions are detailed.

$$ \left(\frac{\partial p}{\partial n}\right)_b = \begin{cases} 0 & \text{if } \dot{d} = 0 \\ -(\rho(\dot{d} \cdot n)) & \text{otherwise} \end{cases}. \qquad (4) $$

<table><caption>**Table 3** Free-stream conditions, instrumented annular ring</caption><tr><td>M</td><td>11</td></tr><tr><td>P[Pa]</td><td>410</td></tr><tr><td>T[K]</td><td>91.2</td></tr><tr><td>u[m/s]</td><td>2160</td></tr><tr><td>AoA[°]</td><td>30.0</td></tr></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 7** Instrumented annular experimental model and pressure sensors locations</figcaption></figure>
## 3.1 Reference code: MISTRAL-CFD

MISTRAL-CFD [24] is developed at R.Tech since 2001. Various gas models can be run assuming either perfect gas, thermochemical equilibrium gas, or thermo-chemical non-equilibrium gas chemistry. While originally developed for high-speed flows, it can handle flows from subsonic to hypersonic flows. Space integration is completed with the Van Leer-AUSM numerical scheme. Second-order precision is achieved with a Monotonic Upstream-centered Scheme for Conservation Laws (MUSCL) extrapolation complemented by the Minmod limiter (Fig. 4). Time integration using a local time step is performed with a Generalized Minimal RESidual (GMRES) method to reach a steady-state solution. The massively multi-blocked structured mesh approach allows flexibility in the geometrical modeling combined with the efficient and accurate resolution of the physical phenomena. The large number of blocks also makes possible the effective domain decomposition for parallelization on large-scale computational facilities. The direct integration of the grid generator GridPro [30] allows the efficient creation of these massively multi-blocked grids.

MISTRAL has been used to create an extensive database for the Debrisk code comprising over 3000 CFD computations [31]. MISTRAL-CFD was initially developed for the aerothermal analysis of spacecraft re-entry vehicles such as IXV and planetary re-entry capsules. These bodies are typically blunt body configurations. Before applying MISTRAL-CFD to the objects available in DEBRISK (boxes, cylinders, rings, sphere caps, etc) for which the flow physics can be significantly more complex than for re-entry vehicles (shock-shock interaction, attached shock, internal shock waves), an extensive validation campaign has been carried out [32].

## 3.2 Presentation of the cases

In this subsection, study cases are briefly presented. Several objectives are pursued to evaluate the BLIZZARD solver: assess the mass flow loss, validate by simulating wind tunnel

<figure><img src="image_10.png" /><figcaption>Fig. 8 Illustration of the attitude of the simulation</figcaption></figure>

<table><caption>**Table 4** Free-stream conditions, Jason satellite</caption><tr><td>M</td><td>26.17</td></tr><tr><td>Altitude[km]</td><td>78</td></tr><tr><td>P[Pa]</td><td>1.1471</td></tr><tr><td>T[K]</td><td>202</td></tr><tr><td>u[m.s<sup>-1</sup>]</td><td>7457.2</td></tr><tr><td>AoA[°]</td><td>135.0</td></tr><tr><td>AoS[°]</td><td>45.0</td></tr></table>

tests, and test an re-entry cases. For each case, the geometry model and mesh are described, and the inflow conditions are given.

### 3.2.1 Sphere

Results concerning the sphere are presented in 4.1. In this paper, the sphere is used for the estimation of the mass loss of the present numerical method, the validations of BLIZZARD on the sphere among other simple geometries can be found in [33]. The model is a 3D sphere of 1 m diameter immersed in a hypersonic flow. The BLIZZARD model is a surface mesh composed of 2814 triangles and is presented Fig. 5. Inflow conditions are detailed in Table 1. BLIZZARD results have been obtained with the steady-state solution considering a perfect gas and an adiabatic wall. For mesh convergence, four wall mesh sizes are used: $[4.5e^{-2}, 2.25e^{-2}, 1.125e^{-2}, 5.625e^{-3}]$ m, giving an approximate diameter ratio $[\frac{D}{22}, \frac{D}{44}, \frac{D}{88}, \frac{D}{176}]$, respectively.

### 3.2.2 VKI static models

Experiments were carried out by the VKI in the Longshot wind tunnel for a test campaign aiming to introduce more complex shapes within the object-oriented destructive re-entry code DEBRISK [32]. During the campaign, both instrumented static runs and free-flight runs were conducted. In the following, comparisons to static cases are presented. Simulations with BLIZZARD of free-flying cases are presented in [20, 21].

For the experiments, the geometries of interest were identified based on the specific flow features they could exhibit for some re-entry attitudes. Attached shock waves, shock-shock interactions, flows over concave surfaces, or unsteady flow configurations were deemed particularly interesting. Among a range of geometries, two were selected for the experimental investigations: a hollow hemisphere (typically representing a propellant reservoir split in two) and an annular ring with an aspect ratio of 1/4th (corresponding to a rocket interstage fairing).

**Instrumented hollow hemisphere** The numerical simulation of the static hemisphere test is shown in sec.5.1. The

<figure><img src="image_11.png" /></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 10** Cube mesh</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Table 5** Initial free-stream conditions and, cube attitude and rotation rates</figcaption></figure>

<table><tr><td>T[K]</td><td>213.78</td></tr><tr><td>P[Pa]</td><td>3.71</td></tr><tr><td>u[m/s]</td><td>4011.77</td></tr><tr><td>Pitch [°]</td><td>66.88</td></tr><tr><td>Yaw [°]</td><td>61.46</td></tr><tr><td>Roll [°]</td><td>0.96</td></tr><tr><td>rotation rate x[°.s<sup>-1</sup>]</td><td>8.3</td></tr><tr><td>rotation rate y[°.s<sup>-1</sup>]</td><td>-2.04</td></tr><tr><td>rotation rate z[°.s<sup>-1</sup>]</td><td>3.4</td></tr></table>

fixed hollow hemisphere is instrumented with 15 pressure sensors, 13 of them located in the concave region. The experimental model and sensor locations are shown in Fig. 6.

Both steady and unsteady time solutions were computed with BLIZZARD. Inflow conditions are detailed in Table 2; note that these are averaged values of the time-dependent real inflow conditions of the experiment. The evolution of the inflow condition was taken into account for BLIZZARD unsteady computations.

**Instrumented annular ring** The numerical simulation of the static ring test is shown in sec. 5.2. The fixed annular ring is instrumented with 15 pressure sensors, 8 of them located in the internal part of the ring, a region where potential shock interactions might occur. The experimental model and sensor locations are shown in Fig. 7.

Averaged experimental free-stream conditions are presented in Table 3. Again, both steady and unsteady solutions was computed with BLIZZARD, and the real evolution
<table><caption>**Table 6** Summary of the test cases</caption><thead><tr><th>Geometry</th><th>Test case type</th><th>Comparison to MISTRAL</th><th>Comparison to the Newtonian theory</th></tr></thead><tbody><tr><td>Sphere</td><td>Mass flow loss assessment</td><td></td><td></td></tr><tr><td>Hollow hemisphere</td><td>Validation with wind tunnel rebuild</td><td>X</td><td></td></tr><tr><td>Annular ring</td><td>Validation with wind tunnel rebuild</td><td>X</td><td></td></tr><tr><td>JASON satellite</td><td>Mass flow loss assessment and application to re-entry case</td><td>X</td><td>X</td></tr><tr><td>Cube</td><td>Rebuild re-entry trajectory</td><td></td><td>X</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Fig. 11** Slice of the BLIZZARD mesh for the coarse and fine mesh, sphere case</figcaption></figure>

of the flow conditions was considered for the unsteady computation.

### 3.2.3 JASON satellite

The BLIZZARD model has been applied on the real shape of the satellite JASON. Jason-1 mission marks the inception of a lasting oceanic observatory, with its measurements benefiting a wide range of earth sciences disciplines, including oceanography, ocean forecasting, marine meteorology, geophysics, and climatology.

This case is particularly noteworthy because, in practical terms, it represents the typical kind of geometry to be simulated with destructive re-entry codes. This type of geometry presents a substantial challenge for traditional CFD codes. The complexity of the geometry becomes a bottleneck for engineers, given that considerable time is required for mesh design, often spanning from several weeks to several months. Furthermore, its unconventional shape introduces complex physical phenomena, such as shock interactions, which are currently not captured by existing spacecraft-oriented codes.

In sec.4.2, the mass loss of the numerical method on this complex shape is studied. Then in sec.6.2, mesh convergence

and comparisons between BLIZZARD, MISTRAL and the Newtonian theory are presented.

Previous MISTRAL-CFD simulations were conducted for a trajectory point located at 78 km in the scope of a contract between R.Tech and CNES for the validation of the PAM- PERO code [22]. The reference only details the comparison with the high-fidelity code Mistral-DSMC [34], a direct simulations Monte Carlo (DSMC) code; however, similar computations in the continuum regime have been carried out by Mistral-CFD.

The velocity is 7457.2 m/s, resulting in a Mach number of 26.17 at this altitude. The Reynolds number for these simulations is $2e^5$, with a Knudsen number of $1.9e^{-4}$ for a reference length of 1 m, indicating laminar flow with slip conditions. The attitude carried out for the simulation is illustrated in Fig.8, with the inflow vector heading into the page. The free-stream conditions are shown in Table 4, where AoA refers to the angle of attack and AoS to the slideslip angle.

The MISTRAL structured multi-block mesh was made using GridPro. It contains more than 14.5 million cells and the height of the first cell layer at the wall is $2e^{-6}$ m (Fig. 9).

For BLIZZARD simulations, the surface mesh input is the MISTRAL body-fitted surface mesh. Four different wall cell sizes were set at the wall at the initialization. From the coarsest to the finest [$5.816e^{-2}$, $2.908e^{-2}$, $1.454e^{-2}$, $7.27e^{-3}$] m.

<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. 12** Surface distribution of dimensionless normal mass flow (top scaled from $-1.5e^{-1}$ to $1.16e^{-1}$, bottom scaled from $-2.84e^{-2}$ to $3.65e^{-2}$), sphere case</figcaption></figure>

### 3.2.4 Cube

In sec.6.1, a comparison between the evolution of the attitude of a cube computed with pressure loads of the modified Newtonian theory and BLIZZARD, is presented. This comparison aims to see the effect of the higher fidelity of the fluid modeling on the evolution of the attitude along a re-entry trajectory. For the simulation of the trajectory, the beginning of the re-entry of a cube has been computed with analytical models, i.e., the modified Newtonian theory for the continuum regime, coupled with an atmosphere model and a 6-DOF rigid body solver. Once the cube is fully in the continuum regime, the same free-stream conditions, which are time-dependent, are imposed as BLIZZARD inflow conditions. The initial inflow conditions and the cube's initial attitude and rotation rate for the BLIZZARD computation are given in Table 5. 50 s of the trajectory has been

simulated using the 6-DOF rigid body motion integration coupled with the Euler solver. The cube is 1 m side centered in (0; 0; 0) with a material density of $ρ_s = 2700kg.m^{-3}$, cube mesh is shown in Fig.10.

### 3.2.5 Summary

Table 6 summarizes the application of the different test cases in this article and whether the results are compared to MIS-TRAL or to the Newtonian theory.
## 4 Mass flow losses of the method

The slip boundary condition imposes the normal velocity to vanish at the wall. As mentioned in [18, 19], the use of immersed boundary methods using ghost cells, as in the present GPFM, makes it difficult to conserve the mass at the wall when a hypersonic flow is considered.

To estimate the mass loss induced by the present GPFM, studies on both a simple sphere and a complex industrial object (JASON satellite) were conducted. In sec.6.2, results on the surface pressure distribution and aerodynamic coefficients are analyzed on the satellite to determine how mass loss can impact the results a complex geometry. For both case, mass losses are observed for 4 different mesh sizes, and convergence is verified by observing at the $L_∞$-error norm and the $L_2$-error norm of the mass flow passing through the fluid interface (denoted $||·||_∞$ and $||·||_2$ respectively). Since slip boundary conditions are set for the velocity, only the normal momentum passing through each surface element of area $ds_i$ is verified (Fig. 11).

$$ ||\epsilon||_\infty = \max \left( \frac{\rho_i \mathbf{u}_i \cdot \mathbf{n}_i \, d \, s_i}{\rho_\infty u_\infty} \right) \\ ||\epsilon||_2 = \sum_i^{n_f} \left( \frac{\rho_i \mathbf{u}_i \cdot \mathbf{n}_i \, d \, s_i}{\rho_\infty u_\infty} \right)^2. \qquad (5) $$

### 4.1 Sphere case

In Fig. 12, the momentum due to the normal velocity $\rho_i \mathbf{u}_i \cdot \mathbf{n}_i$ that passes through each surface element $i$ (mass flow) is shown for the four mesh sizes (top scaled from $-1.5e^{-1}$ to $1.16e^{-1}$, bottom scaled from $-2.84e^{-2}$ to $3.65e^{-2}$). Dimensionless values were obtained by dividing the momentum due to normal velocity by the free-stream momentum $\rho_\infty u_\infty$. Note that the presented values were obtained using the

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

<figure><img src="image_18.png" /><figcaption>Fig. 13 Convergence of mass loss at the wall, sphere case</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 14 Mesh at the wall for BLIZZARD finest case</figcaption></figure>
<table><caption>**Table 7** Mesh sizes, Jason satellite</caption><thead><tr><td>dx at wall [m]</td><td>Initial cell count</td><td>Final cell count</td></tr></thead><tbody><tr><td>5.816e-2</td><td>110 224</td><td>241 697</td></tr><tr><td>2.908e-2</td><td>359 613</td><td>1 037 740</td></tr><tr><td>1.454e-2</td><td>1 266 160</td><td>3 701 970</td></tr><tr><td>7.27e-3</td><td>4 719 280</td><td>18 773 300</td></tr></tbody></table>

interpolation from the fluid mesh onto the sphere surface mesh which also generates small precision losses. In Fig. 13, the $L_\infty$-error norm and the $L_2$-error norm (eq.(5)) convergences are presented with a reference value of 0 given that $\mathbf{u} \cdot \mathbf{n} = 0$ is imposed at the wall. On the surface distribution, the maximum positive losses are observed at the 'slip' region ($\theta \approx 90^\circ$) for all mesh refinement. The maximum negative losses are observed on a radial distribution close to the stagnation point. However, a decrease of the mass loss following the mesh size at the wall is observed. For the coarsest mesh size a maximum mass loss of 12% of inflow momentum is observed, whereas for the finest mesh the maximum loss reached is less than 2%. This was also observed on the $L_\infty$-error norm and $L_2$-error norm plots, which follow a 1st-order and a 2nd-order convergence respectively. Thus, similarly to [18], the present GPFM introduces some mass flow loss at the wall, but this loss decreases with mesh refinement.

## 4.2 Industrial configuration: JASON satellite

BLIZZARD results were obtained with the steady-state solution in perfect gas and adiabatic wall. Table 7 shows for each discretization size the initial cell count, i.e., before any AMR, and the final cell count, i.e. after AMR and convergence. The octree mesh generated at the wall interface for the finest refinement is shown in Fig. 14.

In Fig. 15, the momentum due to normal velocity $\rho_i \mathbf{u}_i \cdot \mathbf{n}_i$ that passes through each surface element $i$ (mass flow) is shown for the four mesh sizes. The color bar is in logarithmic scale from $1e^{-8}$ to $1e^{-4}$. Figure 16 plots the convergence of $L_\infty$-error norm and $L_2$-error norm on the mass flow passing through the surface. To compute the error norm, the normal momentum was multiplied by its corresponding element surface area. In the coarsest case, the refinement level was too coarse to mesh the space between the two central boxes of the satellite, leading to errors in the results in this region, see top left mass flow field in Fig. 15). Thus, the coarsest result is not plotted in this graph for legibility reasons. The $L_\infty$-error norm converges with an order lower than one and the $L_2$-error norm converges following the first order. For this complex case, a much lower convergence order is observed than that obtained in the sphere case, where $L_\infty$ follows the 1st order and $L_2$ the 2nd order; see sec.4.1. When examining

the mass loss distribution over the satellite, higher losses seem to be induced by sharp edges, such as the borders of stages or panels, where attached shocks take place. Regions where shocks interfere with the structure appear to generate additional losses.

## 4.3 Discussions on mass loss

The mass loss generated by the present IBM was studied on two cases, a sphere and a satellite. In both cases, the total mass loss decreased with the mesh refinement. However, the order of convergence differed between the two cases. The mass loss on the sphere follows a 1st-order of convergence for the $L_\infty$ error and 2nd-order for the $L_2$ error whereas that of the satellite followed approximately a 0.5-order and 1st-order respectively. The lower convergence order on the satellite can be explained by the presence of attached shocks and shock interactions on the surface that generate additional mass loss. Sec.6.2 analyses, the comparison made between the present solver and MISTRAL-CFD on the JASON satellite to determine whether the reduction of the mass loss leads to a convergence on the aerodynamic coefficients and on the surface pressure distribution, which are the quantities of interest to compute a re-entry trajectory.

# 5 Simulation of hypersonic wind tunnel static tests

A previous publication [25] studied the validation of the BLIZZARD solver for classical blunt body cases (sphere, diamond airfoil, atmospheric re-entry demonstrator (ARD)). This article aims to validate the BLIZZARD solver for more complex geometries such as hollow hemispheres and annular rings by simulating two wind tunnel experiments conducted by Von Karman Institute for Fluid Dynamics (VKI) [32]. In this section, the experimental data, the results obtained with BLIZZARD and with the body-fitted Navier-Stokes CFD solver MISTRAL-CFD are compared. For the Mistral-CFD results, the mean of the unsteady solution is calculated at pressure sensors, and a snapshot of the unsteady solution is used for the flow and surface fields comparisons. Discrepancies due to viscous effects and the chemistry modeling are expected; these errors have been evaluated by comparing BLIZZARD, MISTRAL-CFD Euler and MISTRAL-CFD Navier-Stokes simulations on various geometries [33, 35]. To summarize, on pressure distribution, less than 5% of error is observed compared to MISTRAL-CFD Euler and around 10% compared to MISTRAL-CFD Navier-Stokes. An error of around 5% was observed on the aerodynamic coefficients. Globally, these discrepancies increase with the
<figure><img src="image_20.png" /><figcaption>Fig. 15 Convergence of mass flow at the wall, JASON case</figcaption></figure>

<figure><img src="image_21.png" /></figure>
<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 16** Mass loss convergence plot, JASON case</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 17** Slice view of the steady solution mesh, X-Y plane (left), X-Z plane (right)</figcaption></figure>

Mach number, and from Mach 20, the viscous effect and the chemistry have a significant impact on the pressure.

## 5.1 Instrumented hollow hemisphere

A slice view of the mesh in the X-Y plane and the X-Z plane of the steady computation is shown in Fig. 17. Due to flow re-circulation in the concavity, the AMR process appears to refine the mesh randomly in this region, however, the stagnation point is still captured at the bottom of the concavity. The mesh of the converged solution contains around 1.2e6 cells (Fig. 17).

Figure 18a shows the front view of BLIZZARD's surface pressure coefficient. In Fig. 18b, surface pressure coefficients obtained at sensor locations along the longitudinal axis with the experiment (red markers), Mistral-CFD (light blue markers), and BLIZZARD (green markers: steady, dark blue:

unsteady) are presented. Fluctuation bars were added for each sensor for both BLIZZARD unsteady computation and the experimental values. These fluctuations are caused by recirculation of the flow in the concave region. Markers with a null pressure coefficient correspond to sensors located at the back face of the model. The experimental measurements fluctuate significantly (30% at sensor located at $α = -70°$), but the mean values remain almost constant in the concave region. Mistral-CFD also computes almost constant mean pressure at all concave-located sensors. The present steady solution appears to slightly underestimate the mean measured pressure at all sensors whereas the unsteady solution is in good agreement with Mistral-CFD results. Snapshots of the BLIZZARD and the experimental Schlierens are shown in Fig. 19. The numerical Schlieren was obtained by computing the gradient of the density field for the unsteady simulation. The numerical view is at $t = 1.5ms$.
<figure><img src="image_25.png" /><figcaption>(b) Wall pressure at sensors along the longitudinal axis. Green crosses: present steady solution, dark blue squares: present unsteady solution, light blue triangles: Mistral-CFD, red dots: experiment</figcaption></figure>

<caption>**Fig. 18** Wall pressure coefficients, instrumented hollow hemisphere case</caption>

which corresponds to the time from which the flow around the hemisphere remains quasi-steady, the first 1.4ms of simulation being characterized by recirculation phenomena in the concave region. The shock position computed at this time is similar to that of the experiment.

## 5.2 Instrumented annular ring

A slice view of the mesh in the X-Y plane and the X-Z plane of the steady computation is shown in Fig. 20. In this case, the front shock is not properly captured with the AMR processing, given that the refinement occurs mainly on the shock interactions with the back face of the ring. In Fig. 21, the near-wall density field of the steady solution is compared with a snapshot of the Mistral-CFD unsteady solution. Similar fields are computed especially at the front face, even

if the present solver does not capture the shock interactions with the back face of the ring with as much precision as the Mistral-CFD solver.

The comparison of all pressure coefficients at the sensors is shown 22. Sensor markers where $\beta < 0$ correspond to sensors located on the front face of the annular ring whereas $\beta > 0$ are sensors located on the back face. At all front sensors present results, either steady or unsteady, are in very good agreement with the experiment and the CFD reference. At sensors located on the back, the present solver does not compute the same values as the experiment or the CFD reference, especially when $\beta > 45$, corresponding to the region where shock interactions occur. Figure 23 compares the pressure distribution on the ring's surface with a snapshot of the unsteady solution from MISTRAL-CFD. The field captured at the back face by MISTRAL is very complex and chaotic, whereas BLIZZARD computes an almost steady symmetric field. This shows the need for high-fidelity CFD solvers such as MISTRAL, that are able to capture precisely these chaotic instabilities. However, the BLIZZARD results are considered satisfactory within the scope of its applications.

# 6 Application

BLIZZARD was tested on two application cases. The first application case studies the effect of a higher fidelity of the fluid modeling on the re-entry trajectory. To do so, a part of the trajectory of a cube was computed with the modified Newtonian theory pressure loads and a 6-DOF solver to compare the resulting evolution of the cube's attitude with the attitude computed using the BLIZZARD pressure loads. Previous comparison of BLIZZARD and MISTRAL-CFD Euler simulations on a cube shape show good agreement on the pressure distribution and the aerodynamic coefficients (see Fig. 24 and Table 8). The second application case consists of the verification of the BLIZZARD mesh convergence on the JASON satellite, followed by the comparison of the pressure loads and aerodynamic coefficients using BLIZZARD, MISTRAL-CFD and the modified Newtonian theory.

## 6.1 Impact of fluid modeling on the stability

50 s of a cube's trajectory, driven by the loads of the modified Newtonian theory and a 6-DOF solver were simulated using pressure loads computed by BLIZZARD. The setup and inflow conditions are presented in the sec.3.2. BLIZZARD simulations start at the conditions corresponding to the 150th second of the trajectory using the analytical models.

<figure><img src="image_26.png" /></figure>
<figure><img src="image_27.png" /><figcaption>**Fig. 19** Schlieren comparison between BLIZZARD and the experiment, instrumented hollow hemisphere</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Fig. 20** Slice view of the steady solution mesh, X-Y plane (left), X-Z plane (right)</figcaption></figure>

### 6.1.1 BLIZZARD convergence

Two fluid mesh sizes; 0.049 vs 0.0245 *m* at the wall and 300 000 cells vs 1 300 000 cells after AMR; and three motion integration time steps, *dt* = [0.05, 0.1, 0.5] *s*, were tested for BLIZZARD computation to verify the convergence. The coarse case with *dt* = 0.1*s* took 2 days to run and the fine case with *dt* = 0.05*s* took 2 weeks. Both ran on a single node of the CNES's cluster TREX with an OpenMP multi-threading on 30 processors. Fig. 25 shows the rotation rates on the three axis and the magnitude. Mesh sizes do not seem to significantly affect the resulting motion of the cube given that with the same time step, the coarse and fine meshes give similar results. Concerning the time step convergence, the rotation rates observed with *dt* = 0.1 *s* and *dt* = 0.05 *s* are similar and their magnitude are about 500°/*s*, whereas the rotation rate with *dt* = 0.5 *s* tends to increase dramatically with a final magnitude of approximately 3000°/*s*, which is an order of magnitude higher.

### 6.1.2 Comparison with the modified Newtonian theory

In Fig. 26, three rotation rates driven by the Newtonian theory are compared to the BLIZZARD results. The BLIZZARD rotation rate reaches much higher values, whereas the Newtonian theory tends to reduce the rotation in the same
<figure><img src="image_29.png" /><figcaption>**Fig. 21** Near-wall density field of present steady solution and MISTRAL-CFD solution</figcaption></figure>

time range. In Fig. 27, snapshots of the first 14 s (every 1 s) of the cube motion computed with the Newtonian theory and BLIZZARD are compared. The BLIZZARD results shown here were computed with the fine mesh and $dt = 0.05$ s. The cube's initially exposed faces can be monitored with the white and the pink markers. With the Newtonian theory, the cube follows the initial rotation rate with a constant motion, whereas with BLIZZARD the cube tends to oscillate with a slight side rotation keeping the white marked face in the front for the first 10s, after which the pink marked face becomes the mainly exposed face. This behavior highlights the capabilities of BLIZZARD to compute "stabilized" motions where the modified Newtonian theory couldn't due to simplified aerodynamic coefficient computation. This observation is very promising and encourages the use of

high-fidelity fluid modeling to properly assess debris trajectories and demise.

## 6.2 Surface loads and aerodynamic coefficient estimation on an industrial case

The setup and inflow conditions are presented in sec.3.2.

### 6.2.1 BLIZZARD mesh convergence

BLIZZARD results were obtained with the steady-state solution considering a perfect gas and an adiabatic wall. The drag coefficient $C_D$ and the max surface pressure coefficient $C_{pmax}$ are observed to determine the convergence of the cases. For both $C_D$ and $C_{pmax}$, a variation under 0.2% was set to trigger the end of the AMR process in BLIZZARD.

<figure><img src="image_30.png" /></figure>
<figure><img src="image_31.png" /><figcaption>**Fig. 22** Wall pressure coefficient at sensors. Green crosses: present steady solution, dark blue squares: present unsteady solution, light blue triangles: MISTRAL-CFD, red dots: experiment</figcaption></figure>

**Surface pressure:** Figure 28 shows the mesh convergence on the distribution of the surface pressure coefficient. The camera is aligned with the flow direction and the windward side is visible. For all coarsening, the results compute the interactions between the shock generated by the solar panel and the upper stage and the front parabola. The more the mesh is refined, the more accurately is the additional pressure caused by shock interactions captured.

**Integrated coefficients:** In Table 9, drag, lateral lift, lift, and pitching moment coefficients are presented for each mesh refinement. The reference length used is $L_{ref} = 1m$ and the reference surface is $S_{ref} = 0.327m^2$ corresponding to the projected area of the satellite at the computed attitude ($AoA = 135°$, $AoS = 45°$). Coefficients are computed as follows:

<figure><img src="image_32.png" /><figcaption>**Fig. 23** BLIZZARD (steady solution) and MISTRAL-CFD (snapshot of unsteady solution) surface pressure coefficients distributions, instrumented annular ring</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>(a) Surface pressure coefficient on the exposed face, comparison between BLIZZARD (top) and MISTRAL-CFD Euler (bottom)</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>(b) $C_p^{max}$ distribution over the curve cutting the cube in y direction. Black: modified Newtonian theory, green: BLIZZARD, blue: MISTRAL-CFD Euler</figcaption></figure>

<caption>**Fig. 24** Comparison of BLIZZARD and MISTRAL-CFD Euler on the pressure distribution on a cube</caption>
<table><caption>**Table 8** Comparison of drag coefficient and $C_p^{max}$, cube shape</caption><thead><tr><th>Solver</th><th>C<sub>D</sub></th><th>C<sub>p</sub><sup>max</sup></th></tr></thead><tbody><tr><td>MISTRAL-CFD Euler</td><td>1.643</td><td>1.818</td></tr><tr><td>BLIZZARD</td><td>1.712</td><td>1.823</td></tr><tr><td>Newtonian theory</td><td>1.77</td><td>1.821</td></tr></tbody></table>

$$C_D = \frac{F_x^i}{\rho_\infty ||\mathbf{u}_\infty|| S_{ref}}, \qquad (6)$$

$$C_L^{lat} = \frac{F_y^i}{\rho_\infty ||\mathbf{u}_\infty|| S_{ref}}, \qquad (7)$$

$$C_L = \frac{F_z^i}{\rho_\infty ||\mathbf{u}_\infty|| S_{ref}}, \qquad (8)$$

$$C_m = \frac{M_z^b}{\rho_\infty ||\mathbf{u}_\infty|| l_{ref} S_{ref}}, \qquad (9)$$

here $F_x^i, F_y^i, F_z^i$ are the aerodynamic integrated force components in the aerodynamic reference frame and $M_z^b$ is the integration pitch moment in the solid reference frame. In Fig. 29, the convergence relative to the finest mesh size for the four coefficients is plotted. The drag plot shows that the coarsest mesh results in a closer $C_D$ than the medium-coarse ($dx = 2.91e-2$ m) and medium-fine ($dx = 1.45e-2$ m) cases. This is due to errors in the results in the region between the central boxes of JASON, where the refinement level of the coarsest mesh has a cell size too large to mesh the space between the boxes. This has a lower impact on the other aerodynamic coefficients, where a convergence in accordance with the mesh size is observed, following a order between 1 and 2.

### 6.2.2 Comparison with MISTRAL-CFD and Newtonian theory

The MISTRAL-CFD simulations were performed solving the Navier-Stokes equations under chemical and thermal non-equilibrium conditions, with the laminar solver and with a wall temperature set at 1200K.

The durations of input mesh setup and of running the simulations are 2 month of work to create the structured multi-block mesh and 1 month of computation of 200 CPU (144 000 cpu hours) to obtain the MISTRAL results. With the modified Newtonian theory, the mesh is easier to design and was done in 1 day and the surface pressure distribution and aerodynamical coefficients being computed within seconds. For BLIZZARD, the same meshing strategy as

spacecraft-oriented code was implemented, which required 1 day of engineering work. Runs were launched on 1 CPU with 32 OpenMP multi-threading; for the coarsest mesh, the steady solution was obtained within a couple of hours while it took 2 days for the finest mesh to converge.

Figure 30 displays the Mach field and the corresponding numerical Schlieren obtained for the finest BLIZZARD case. Shocks are accurately captured near the object. What's noteworthy is the intricate flow patterns generated, where interactions between shocks and structures are evident, which means that the shocks originating from the panels interact with those emanating from the structure at the rear of the satellite, resulting in numerous shock-shock interactions.

Figure 31 provides visual representations of the surface pressure coefficient $C_p$ distribution across the object's surface. BLIZZARD results are compared with those computed through MISTRAL-CFD and a spacecraft-oriented code. For this attitude, there are multiple shock-shock interactions which results in a complex flow structure. Obtaining comparable results with BLIZZARD presents a challenge, given that there are discrepancies with the reference regarding the locations of the shock interactions. This case represents a configuration where spacecraft-oriented codes produce significantly inferior outcomes, especially on the upper stage box, whereas BLIZZARD captures a similar interaction to that of MISTRAL-CFD.

In Table 10, the aerodynamic coefficients as well as the $C_p^{max}$ obtained by the three codes are compared. With the present solver, a 1.19% difference for $C_D$ is obtained, compared to 23.82% with the spacecraft-oriented code. The discrepancies on the $C_L$ are 5.26% compared to 19.6%. For the $C_p^{max}$, BLIZZARD value is lower than MISTRAL one, but still significantly improved compared to the modified Newtonian theory.

For a better assessment of the impact on the aerodynamic coefficients and the wall pressure when using BLIZZARD, two other orientations of the satellite are compared. Additional simulations were run at $AoA = 0°$, $AoS = 0°$ (Table 11) $AoA = -70°$, $AoS = 0°$ (Table 12). For $AoA = 0°$, $AoS = 0°$, a discrepancy with reference of 4.2% is observed on the drag compared to 30% the Newtonian theory, 5.3% on the lift compared to 25%. For $AoA = -70°$, $AoS = 0°$, differences of 1.3% are computed on the $C_D$ versus 10% with the analytical method, 7.3% versus 13.9% on the $C_L$. For both cases, BLIZZARD underestimates the maximum wall pressure compared to MISTRAL, but significant improvements are still observed compared to the Newtonian theory.
<figure><img src="image_35.png" /><figcaption>**Fig. 25** BLIZZARD cube rotation rate convergence</figcaption></figure>
<figure><img src="image_36.png" /><figcaption>**Fig. 26** Cube rotation rate evolution with modified Newtonian theory loads</figcaption></figure>

# 7 Conclusion

The re-entry CFD solver BLIZZARD aims to bridge the gap between analytic/engineering methods (spacecraft-oriented methods) and 'high-fidelity' methods (high-resolution body-fitted CFD) for pressure load assessment in the scope of the implementation of the high-fidelity branch of the spacecraft-oriented tool PAMPERO. BLIZZARD solves the Euler equations closed by the perfect gas relation. The fluid domain meshing is fully automated using an Octree method and it is coupled to a ghost-cell IBM for the near-wall modeling. AMR methods are implemented to optimize the cell count while refining the regions of interest (shock, wake, shock interactions).

The mass loss generated by a ghost-cell IBM is a classic drawback of the IBM method [18]. In the simulations by BLIZZARD on a sphere and on the JASON satellite a reduction of the mass loss consistent with the mesh size is observed. However, the convergence order is higher in the sphere case than in the satellite case given that the satellite's attitude is such that at some locations, the shocks are

attached to the surface and interact with the geometry, creating additional mass loss.

To validate BLIZZARD, simulations of instrumented VKI wind tunnel static tests were performed, which involve a fixed hollow hemisphere where the concave face is exposed to the flow and a fixed annular ring with an angle of attack of 30°. Concerning the hollow hemisphere, after a few milliseconds of a completely unsteady flow behavior due to recirculation in the concave region, the flow stabilizes and pressure at sensors is in good agreement with both the experiment and the reference solver MISTRAL-CFD. Concerning the ring, no steady solution is obtained due to shock interactions with the back face. With the reference CFD, a strongly chaotic behavior of the flow on the back face is observed, whereas BLIZZARD computes an almost symmetric field, probably due to the cell size being too large to accurately capture this complex flow. At pressure sensors located on the front face, the experiment, MISTRAL-CFD and BLIZZARD are in good agreement, whereas for the sensors located at the back face, more discrepancies are observed.

BLIZZARD was tested on two re-entry cases, a cube and the JASON satellite. Simulations of part of a cube's

<figure><img src="image_37.png" /></figure>
<figure><img src="image_38.png" /><figcaption>Fig. 27 Newtonian theory vs BLIZZARD cube motion 0 to 14 s</figcaption></figure>
<figure><img src="image_39.png" /><figcaption>**Table 9** Drag, lateral lift, lift and pitching moment coefficients for each mesh size, Jason satellite</figcaption></figure>

<table><thead><tr><th>dx at wall [m]</th><th>C<sub>D</sub></th><th>C<sub>L</sub><sup>lat</sup></th><th>C<sub>L</sub></th><th>C<sub>m</sub></th></tr></thead><tbody><tr><td>5.816e-2</td><td>1.6545</td><td>0.4848</td><td>0.6806</td><td>0.00321</td></tr><tr><td>2.908e-2</td><td>1.6392</td><td>0.4612</td><td>0.8249</td><td>0.01454</td></tr><tr><td>1.454e-2</td><td>1.6494</td><td>0.4346</td><td>0.8433</td><td>0.01245</td></tr><tr><td>7.27e-3</td><td>1.6605</td><td>0.4273</td><td>0.8631</td><td>0.01260</td></tr></tbody></table>

trajectory initially computed using pressure loads from the modified Newtonian theory were performed by BLIZZARD. The induced rotation rates are 2 orders of magnitude higher than those computed with the analytical method. In addition, a better pressure distribution computed with a higher-fidelity modeling allows a stabilization of the cube's attitude that cannot be obtained with the modified Newtonian theory.

The pressure loads and aerodynamic coefficients computed by BLIZZARD, MISTRAL-CFD, and the modified

Newton theory on the JASON satellite were compared. After the verification of BLIZZARD's mesh convergence, a significant improvement on the computed surface pressure field was observed compared to that of the modified Newtonian theory. Similar shock interactions with the surface are captured by BLIZZARD and MISTRAL-CFD. The aerodynamic coefficients computed by BLIZZARD are in good agreement with those computed by MISTRAL-CFD, but not with those of the modified Newtonian theory (1–5% errors for BLIZZARD and ≈ 20% errors for the analytical method). In addition to the improved results compared to the modified Newtonian theory, BLIZZARD needs much less setup and computation time than MISTRAL-CFD (3–4 months for MISTRAL-CFD and 4 days for BLIZZARD's finest case).

This article highlighted the improvements brought by BLIZZARD in the pressure load and aerodynamic coefficient estimation compared to the analytical formulation currently used in existing spacecraft-oriented codes. New methods are being developed for BLIZZARD in order for it to be
<figure><img src="image_40.png" /><figcaption>Fig. 29 Integrated coefficients convergence relative to the finest case</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>Fig. 30 Mach distribution (left) and corresponding schlieren (right)</figcaption></figure>
<figure><img src="image_42.png" /><figcaption>**Fig. 31** Comparison of the surface pressure coefficient distributions</figcaption></figure>

integrated in the high-fidelity version of the re-entry tool PAMPERO and used for a full re-entry application. These methods include fluid-structure coupling with a FEM solver to compute the local stress distribution and structure failure using high-fidelity methods, and Machine Learning modules to assess surface convective heat fluxes and viscous effects.

<table><caption>**Table 10** Comparison of the aerodynamic coefficient and $C_p^{max}$, Jason satellite, AoA = 135°, AoS = 45°</caption><thead><tr><th>Solver</th><th>C<sub>D</sub></th><th>C<sub>L</sub><sup>lat</sup></th><th>C<sub>L</sub></th><th>C<sub>p</sub><sup>max</sup></th></tr></thead><tbody><tr><td>MISTRAL-CFD</td><td>1.641</td><td>0.387</td><td>0.820</td><td>8.99</td></tr><tr><td>BLIZZARD</td><td>1.6605</td><td>0.4273</td><td>0.8631</td><td>5.52</td></tr><tr><td>Newtonian theory</td><td>1.212</td><td>0.405</td><td>0.615</td><td>1.831</td></tr></tbody></table>

<table><caption>**Table 11** Comparison of the aerodynamic coefficient and $C_p^{max}$, Jason satellite, AoA = 0°, AoS = 0°</caption><thead><tr><th>Solver</th><th>C<sub>D</sub></th><th>C<sub>L</sub><sup>lat</sup></th><th>C<sub>L</sub></th><th>C<sub>p</sub><sup>max</sup></th></tr></thead><tbody><tr><td>MISTRAL-CFD</td><td>1.742</td><td>2.22e-3</td><td>0.102</td><td>8.20</td></tr><tr><td>BLIZZARD</td><td>1.693</td><td>4.38e-3</td><td>0.098</td><td>3.74</td></tr><tr><td>Newtonian theory</td><td>1.425</td><td>1.27e-5</td><td>0.123</td><td>1.831</td></tr></tbody></table>

Works have also been initiated to simulate the rarefied and transitional regime with BLIZZARD.
<table><caption>**Table 12** Comparison of the aerodynamic coefficient and $C_p^{max}$, Jason satellite, AoA = -70°, AoS = 0°</caption><thead><tr><th>Solver</th><th>C<sub>D</sub></th><th>C<sub>L</sub><sup>lat</sup></th><th>C<sub>L</sub></th><th>C<sub>p</sub><sup>max</sup></th></tr></thead><tbody><tr><td>MISTRAL-CFD</td><td>4.450</td><td>6.305e-4</td><td>1.144</td><td>22.56</td></tr><tr><td>BLIZZARD</td><td>4.392</td><td>3.014e-4</td><td>1.227</td><td>14.07</td></tr><tr><td>Newtonian theory</td><td>3.971</td><td>1.096e-4</td><td>1.303</td><td>1.831</td></tr></tbody></table>

**Author Contributions** V.L. and E.C wrote the main manuscript. I.M. contributed to BLIZZARD model research and development and also to some of BLIZZARD's results P.V.H. and I.M. produced MISTRAL results. All authors reviewed the manuscript.

**Funding** Open access funding provided by Centrale Nantes.

**Data Availability** No datasets were generated or analyzed during the current study.

## Declarations

**Conflict of interest** The authors declare no conflict of interest.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

1. De Persis, C., Lemmens, S.: An overview of the destructive re-entry analysis tools available in europe: State of the art, latest advances, and open points. J. Space Saf. Eng. **11**(1), 35–43 (2024). https://doi.org/10.1016/j.jsse.2023.11.009

2. Peddakotla, S.A., Morgado, F., Thillaithevan, D., O Driscoll, D., Santer, M., Maddock, C., Vasile, M., Fossati, M.: A Multi-fidelity and Multi-disciplinary Approach for the Accurate Simulation of Atmospheric Re-entry. 73th International Astronautical Congress (2022)

3. Bariselli, F.: Development of an Enhanced Spacecraft Fragmentation Code. ATD3 Workshop (2024)

4. Peskin, C.S.: Flow patterns around heart valves: A numerical method. J. Comput. Phys. **10**(2), 252–271 (1972). https://doi.org/10.1016/0021-9991(72)90065-4

5. Roma, A.M., Peskin, C.S., Berger, M.J.: An adaptive version of the immersed boundary method. J. Comput. Phys. **153**(2), 509–534 (1999). https://doi.org/10.1006/jcph.1999.6293

6. Pinelli, A., Naqavi, I.Z., Piomelli, U., Favier, J.: Immersed-boundary methods for general finite-difference and finite-volume navier stokes solvers. J. Comput. Phys. **229**(24), 9073–9091 (2010). https://doi.org/10.1016/j.jcp.2010.08.021

7. Constant, E., Favier, J., Meldi, M., Meliga, P., Serre, E.: An immersed boundary method in openfoam: verification and validation. Comput. Fluids **157**(3), 55–72 (2017)

8. Hamza, R., Constant, E., Favier, J., Meliga, P., Serre, E., Meldi, M., Goncalves, E.: Direct Numerical Simulation of Compressible Flows Around Spherical Bodies Using the Immersed Boundary Method, pp. 325–331 (2019). https://doi.org/10.1007/978-3-030-04915-7_43

9. Kumar, V., Sharma, R.K.S.A.: Central upwind scheme based immersed boundary method for compressible flows around complex geometries. Computers and Fluids (2020)

10. Kim, J., Kim, H.C.D.: An immersed-boundary finite-volume method for simulations of flow in complex geometries. Journal of Computational Physics (2001)

11. Ren, W., Shu, W.Y.C.: An efficient immersed boundary method for thermal flow problems with heat flux boundary conditions. International Journal of Heat and Mass Transfer (2013)

12. Park, S.G., Chang, C.B., B., K., Sung, H.J.: Simulation of fluid flexible body interaction with heat transfer. International Journal of Heat and Mass Transfer (2017)

13. Schneiders, L., Günther, C., Meinke, M.: An efficient conservative cut-cell method for rigid bodies interacting with viscous compressible flows. J. Comput. Phys. **311**, 62–86 (2016). https://doi.org/10.1016/j.jcp.2016.01.026

14. Muralidharan, B., Menon, S.: Simulation of moving boundaries interacting with compressible reacting flows using a second-order adaptive cartesian cut-cell method. J. Comput. Phys. **357**, 230–262 (2018). https://doi.org/10.1016/j.jcp.2017.12.030

15. Cui, Z., Yang, Z., Jiang, H.-Z.: A sharp-interface immersed boundary method for simulating incompressible flows with arbitrarily deforming smooth boundaries. Int. J. Comput. Methods **15**(01), 1750080 (2017). https://doi.org/10.1142/s0219876217500803

16. Qu, Y., Shi, R., Batra, R.-C.: An immersed boundary formulation for simulating high-speed compressible viscous flows with moving solids. J. Comput. Phys. **354**, 672–691 (2018). https://doi.org/10.1016/j.jcp.2017.10.045

17. Khalili, M.E., Larsson, M., Müller, B.: Immersed boundary method for viscous compressible flows around moving bodies. Comput. Fluids **170**, 77–92 (2018). https://doi.org/10.1016/j.compfluid.2018.04.033

18. Baskaya, A.O., Capriati, M., Turchi, A., Magin, T., Hickel, S.: Assessment of immersed boundary methods for hypersonic flows with gas surface interactions. Comput. Fluids **270**, 106134 (2024). https://doi.org/10.1016/j.compfluid.2023.106134

19. Fadlun, E.A., Verzicco, R., Orlandi, P., Mohd-Yusof, J.: Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations. J. Comput. Phys. **161**(1), 35–60 (2000). https://doi.org/10.1006/jcph.2000.6484

20. Ledermann, V.: Modelisation fluide-structure d'une rentr e atmospherique d'un d bris spatial. PhD thesis, Central Nantes (2025)

21. Ledermann, V., Constant, E., Spel, M., Stainier, L., Annaloro, J., Galera, S.: Pampero hifi : improvements in mechanical and fluid-structure interaction modelling for a spacecraft-oriented code. CEAS Space Journal (2025)

22. Hauwaert, P.V., Constant, E., I., M., Spel, M., Makowski, P., Ledermann, V., Dumon, J., Galera, S., Annaloro, J.: Pampero v3, a spacecraft-oriented reentry analysis code. In: 2nd International Conference on Flight Vehicles, Aerothermodynamics and Re-entry Missions Engineering (FAR) (2022)

23. Grossir, G., Puorto, D., Ilich, Z., Paris, S., Chazot, O., Rumeau, S., Spel, M., Annaloro, J.: Aerodynamic characterization of space debris in the vki longshot hypersonic tunnel using a free-flight measurement technique. Experiments in Fluids **61**(7) (2020). https://doi.org/10.1007/s00348-020-02995-7
24. Scott, T., Dieudonn, W., Spel, M.: Mistral: Conus debut flow field and heat transfer calculations. 34th AIAA Fluid Dynamics Conference and Exhibit (2012)

25. Midani, I.: Maillage automatique et m thodes machine learning pour la pr diction de grandeurs fluides pour les rentr es atmosphriques. PhD thesis, Aix-Marseille Universit (2024)

26. Optimad Engineering Srl: "Bitpit : open source library for scientific high performance computing". Optimad. https://optimad.github.io/bitpit/ (2024)

27. Áfra, A.T., Wald, I., Benthin, C., Woop, S.: Embree ray tracing kernels: Overview and new features. In: ACM SIGGRAPH 2016 Talks. SIGGRAPH '16, pp. 52–1522. ACM, New York, NY, USA (2016). https://doi.org/10.1145/2897839.2927450

28. Péron, S., Benoit, C., Renaud, T., Mary, I.: An immersed boundary method on cartesian adaptive grids for the simulation of compressible flows around arbitrary geometries. Engineering with Computers 37(3), 2419–2437 (2021)

29. Franke, R.: Scattered data interpolation: Tests of some method. Math. Comput. 38(157), 181 (1982). https://doi.org/10.2307/2007474

30. GridPro, Program Development Corporation. http://www.gridpro.com

31. Van Hauwaert, P., Spel, M., Stephane, G., Annaloro, J., Omaly, P.: Cfd based interpolation of drag and heat rates of different objects re- entering the earth atmosphere in a prescribed motion. In: 2nd International Conference on High-Speed Vehicle Science Technology (HISST) (2022)

32. Annaloro, J., Galera, S., Thiebaut, C., Spel, M., Van Hauwaert, P., Grossir, G., Paris, S., Chazot, O., Omaly, P.: Aerothermodynamics modelling of complex shapes in the debrisk atmospheric reentry tool: Methodology and validation. Acta Astronaut. 171, 388–402 (2020). https://doi.org/10.1016/j.actaastro.2020.03.006

33. Midani, I., Ledermann, V., Constant, E., Spel, M., Stainier, L., Boivin, P., Favier, J.: Immersed boundary formulation for complex geometries in hypersonic flows: Application to atmospheric reentry. Comput. Fluids 299, 106715 (2025). https://doi.org/10.1016/j.compfluid.2025.106715

34. Van Hauwaert, P.: Approximation models for drag and heat flux of random tumbling objects in the transitional regime. CEAS Space J. 15(6), 895–921 (2023). https://doi.org/10.1007/s12567-023-00492-1

35. Midani, I., Ledermann, V., Constant, E., Spel, M., Stainier, L., Annaloro, J., Galera, S., Boivin, P., Favier, J., Henneton, M.: Development of a high-fidelity spacecraft-oriented tool. In: 2nd International Conference on Flight Vehicles, Aerothermodynamics and Re-entry Missions Engineering (FAR) (2022)

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.