# Immersed boundary formulation for complex geometries in hypersonic flows: Application to atmospheric reentry

Iko Midani, Valentin Ledermann, Eddy Constant, Martin Spel, Laurent Stainier, Pierre Boivin, Julien Favier

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

## ► To cite this version:

Iko Midani, Valentin Ledermann, Eddy Constant, Martin Spel, Laurent Stainier, et al., Immersed boundary formulation for complex geometries in hypersonic flows: Application to atmospheric reentry. Computers and Fluids, 2025, 299, 10.1016/j.compfluid.2025.106715. hal-05468238

HAL Id: hal-05468238

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

Submitted on 20 Jan 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_3.png" /></figure>

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

# Immersed boundary formulation for complex geometries in hypersonic flows: Application to atmospheric reentry

Iko Midani<sup>a,c,*</sup>, Valentin Ledermann<sup>b,c ID</sup>, Eddy Constant<sup>c</sup>, Martin Spel<sup>c</sup>, Laurent Stainier<sup>b ID</sup>, Pierre Boivin<sup>a ID</sup>, Julien Favier<sup>a</sup>

## ARTICLE INFO

**Keywords:**

Hypersonic flows

Atmospheric re-entry

Immersed boundary methods

Euler equations

## ABSTRACT

This article discusses the challenges of modeling atmospheric reentry using computational fluid dynamics (CFD) due to its complexity and practical industrial applications. Ablation phenomena caused by high energy make it difficult and time-consuming to use a “high-fidelity” CFD method to accurately measure forces and heat flux. As a result, methods based on Newton’s theory are used to model aerodynamic forces, incorporating statistical correlations from CFD results to estimate heat fluxes. However, these methods sacrifice accuracy for CPU and engineering time and have difficulty representing realistic physics when complex phenomena occur, such as shock interactions. To bridge the gap between approximate and high-fidelity methods, we propose a new approach using an automatic grid generation method of the octree Cartesian type coupled to a solver solving Euler’s equations. To apply the boundary condition accurately we compared two immersed boundary methods: a diffuse interface method and a sharp interface method under hypersonic flow configurations. We present a comparative study of these two formulations on verification and validation phases, including an academic test case and an industrial case on a real reentry spacecraft. The novelty lies in applying these methods to complex cases involving strong discontinuities (attached shocks). After concluding this comparative study, we demonstrate that with adapted formulations and an optimized approach, IB methods can handle complex geometries typical of atmospheric reentry.

## 1. Introduction

For many years, modeling atmospheric reentry 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 energy involved in re-entry, and using a “high-fidelity” CFD method to accurately measure forces and heat flux is difficult and time-consuming as it requires manual mesh engineering. As a result, methods based on Newton’s theory are used to model aerodynamic forces, incorporating statistical correlations from CFD results to estimate heat fluxes, while sacrificing accuracy for CPU and engineering time. Moreover, these methods, although extremely efficient in computational time, have much difficulty in representing realistic physics when complex phenomena due to the nature of the geometry as shock interactions appear. To bridge the gap between approximate and high-fidelity methods, manufacturers are increasingly demanding precise modeling of aerothermodynamic phenomena throughout a trajectory with realistic time and engineering costs.

To overcome these issues we developed a new approach using an automatic grid generation method of the octree Cartesian type. The code, called Blizzard uses a solver solving Euler’s equations. For several years, a growing scientific enthusiasm has developed for this type of method which makes it possible to perform CFD simulations on geometries of any kind, where the use of body-fitted meshes is hard to set up. The major drawback of this type of approach comes from the quality of the near-wall modeling due to the nature of Cartesian meshes. Many methods have emerged since to overcome these problems. Among them are the popular, cut-cell, overset, and Immersed Boundary Methods (IBM) chosen for the present work. IBM methods are popular approaches, used initially to model fluid flows around curved and/or moving bodies. These methods were introduced by Peskin, and have been then applied to solve fluid flow problems for many applications [1–4]. Several approaches of IBMs have since been investigated in order to comply with specific applications. In the class

<footnote><sup>a</sup> Aix Marseille Univ, CNRS, Centrale Méditerranée, M2P2, Marseille, France</footnote>

<footnote><sup>b</sup> Nantes Université, CNRS, Ecole Centrale Nantes, GeM, Nantes, France</footnote>

<footnote><sup>c</sup> RTech, Verniolle, France</footnote>

<footnote>* Corresponding author at: Aix Marseille Univ, CNRS, Centrale Méditerranée, M2P2, Marseille, France.
E-mail address: iko.midani@rtech.fr (I. Midani).</footnote>
of direct forcing methods, these approaches could be classified into two categories: the diffuse interface method and the sharp interface method.

In the diffuse-interface IBM the boundary is spread by distributing forcing and source terms [5,6] to the surrounding Cartesian grid cells using a window function. This approach is simpler when it comes to Dirichlet boundary conditions, but it is more difficult to impose Neumann boundary conditions with the required precision. Therefore, the diffuse interface is usually implemented in incompressible flow solvers, in which the velocity components are usually attributed to the location of the boundary. One can find some applications of this method for compressible solvers in the following references: [6–10]. A notable advantage of the diffuse interface IBM is its independence from spatial discretization, making it easier to integrate into existing fluid solvers.

Alternatively, the discrete forcing approaches for tracking boundaries in fluid simulations include the sharp-interface IBM. Several methods fall under this category, such as the cut-cell method [11,12], Cartesian IB method [13], and ghost-cell method [14,15]. 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 IB method, points within the solid body that have at least one neighboring point in the fluid are known as IB 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 IB method but with ghost points representing solid points adjacent to the fluid domain. The reconstruction procedure occurs at these ghost points. While the Cartesian IB method and the ghost-cell method offer less accuracy 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. Since we want to apply the method to complex geometries the ghost cell method has been selected for the study.

In this paper, we conduct a comparative study between these two approaches in the framework of hypersonic flow configurations. In the first part, the governing Euler equations are recalled and we quickly present the associated numerical methods used for their resolution. Secondly, the numerical methods developed for both sharp and diffuse interface methods are exhibited. As verification and validation phases, we propose a complete comparison of these two approaches on several test cases. First, the simulation of a hypersonic flow past a sphere is studied. The application of the two methods is then compared on a real case of an atmospheric reentry vehicle of the blunt body type. Finally, we study the behavior of the methods on a more complex case that includes discontinuities and attached shocks: the diamond airfoil.

By analyzing and comparing these two approaches, we are able to demonstrate that the application of IBM in the context of solving hypersonic flows for complex geometries (including attached shocks) is robust. This point represents a significant advancement as it contributes to demonstrating the interest of this type of approach for industrial applications.

## 2. Governing equations and numerical methods

Euler equations are given as:

$$ \frac{\partial \mathbf{Q}}{\partial t} + \nabla \cdot \mathcal{F}(\mathbf{Q}) = 0 \qquad (2.1) $$

with $\mathbf{Q}$ the conservative variables vector and $\mathcal{F}$ the inviscid flux redtensör:

$$ \mathbf{Q} = \begin{bmatrix} \rho \\ \rho u \\ \rho v \\ \rho w \\ \rho E \end{bmatrix}, \quad \mathcal{F}(\mathbf{Q}) = \begin{bmatrix} \rho u & \rho v & \rho w \\ \rho u^2 + p & \rho u v & \rho u w \\ \rho v u & \rho v^2 + p & \rho v w \\ \rho w u & \rho w v & \rho w^2 + p \\ u(\rho E + p) & v(\rho E + p) & w(\rho E + p) \end{bmatrix} \qquad (2.2) $$

<figure><img src="image_5.png" /><figcaption>Fig. 1. Illustration of the adopted flux convention.</figcaption></figure>

<table><caption>Table 1
Initial conditions on both left and right side of the tube.</caption><thead><tr><th></th><th>Left</th><th>Right</th></tr></thead><tbody><tr><td>p</td><td>10 000 Pa</td><td>1000 Pa</td></tr><tr><td>T</td><td>348 K</td><td>278 K</td></tr><tr><td>u</td><td>0 m s</td><td>0 m s</td></tr></tbody></table>

with, $p$ the pressure, $\rho$ the volume mass, and $(u, v, w)$ the three components of the velocity vector $\mathbf{u}$. $E$ is the total energy and can be written as $E = e + \frac{1}{2}\mathbf{u}^2$ with $e$ is the specific internal energy. The system (2.1) is closed using the ideal gas equation of state:

$$ \begin{cases} p = \rho r T \\ e = c_v T \end{cases} \qquad (2.3) $$

where $c_v$ is the specific heat at constant volume, $r$ the specific gas constant and $T$ the temperature.

The three-dimensional equations (2.1) are solved on a non-uniform Cartesian mesh. For the numerical resolution, we apply finite volume method where the fluid variables are stored at cell centers. The conversion of the volume integration to surface integration is carried out by using the divergence theorem involving the evaluation of fluxes at the faces of each cell. We solve the steady state of (2.1) where the inviscid flux is discretized using a first-order flux vector splitting scheme Van Leer [16] which offers a good compromise between stability and computational cost.

In three-dimension, the spatial term in (2.1) is written in discretized form using the divergence theorem:

$$ \nabla \cdot \mathcal{F}_i(\mathbf{Q}) = \frac{1}{\nu_i} \sum_{j=1}^{n=6} \sum_{k=1}^{m} f_{jk}(\mathbf{Q}) \cdot s_{jk} \qquad (2.4) $$

with $f$ the flux evaluated at the face, the suffix $i$ denote the cell index, $j$ the face index, $k$ the subface index, $\nu$ the volume of the cell, and $s$ the area of each subface. An illustration of the flux convention is shown in Fig. 1.

Temporal integration of (2.1) is carried out by using a first-order explicit upwind scheme. We can express the full discretization in space and time where the inviscid flux is evaluated at time step $n$:

$$ \frac{\mathbf{Q}_i^{n+1} - \mathbf{Q}_i^n}{\Delta t} = -\frac{1}{\nu_i} \sum_{j=1}^{n=6} \sum_{k=1}^{m} f_{jk}(\mathbf{Q}_i^n) \cdot s_{jk} \qquad (2.5) $$

As validation of the solver, an academic test is presented: the Sod shock tube problem [16,17]. The tube is partitioned into two different parts with the initial conditions as follows (see Table 1):

At time $t_0 = 0$ s, the partition is removed and we study the evolution of the solution. Fig. 2 shows the density profile at time $t = 0.6$ with 176 points for the spatial discretization. The computation is performed using a CFL number of 0.8.
<figure><img src="image_6.png" /><figcaption>Fig. 2. Comparison of numerical and analytical density profile at time t = 0.06 s, for the classical SOD configuration [17], red line: numerical, dash line: analytical solution. The numerical solution is obtained for 176 grid points.</figcaption></figure>

## 3. Two immersed boundary methods

In this section, we will develop the methodology for both sharp and diffuse interface IB methods.

**Method 1, Sharp interface:** Imposition of the boundary condition by imposing the primitive variable directly on the boundary cells.

**Method 2, Diffuse interface:** Imposition of the boundary condition by applying a forcing term on the Euler system of equations.

Both use the same classical data structure [18], reminded in Appendix A

### 3.1. Method 1: Sharp interface

In ghost-cell-based sharp interface approaches, the computational domain is divided into fluid, solid, and ghost cells as presented in Appendix A. The values of variables in ghost cells are determined by interpolating from nearby fluid cells and considering the boundary conditions, which impact the flow field through the ghost cells. Ghost cells are solid cells with at least one neighboring fluid cell in any direction (south/north, east/west, or top/bottom). Initially, a simple linear extrapolation method estimates flow variables at ghost cells, but this often results in large and negative weighting coefficients, causing numerical instability with the flow solver.

To address this issue, a more effective procedure introduces a mirror point within the flow domain for each ghost cell, located at some distance from the boundary perpendicular to it. Variables at these mirror points are interpolated from surrounding cells. The interpolation at the image point is crucial, affecting the accuracy and stability of the numerical solution. Direct linear interpolation [18–20] (bilinear in 2D and trilinear in 3D) for second-order reconstruction can cause stability issues near the immersed boundary or when interpolation stencils are insufficient. On the other hand, High-order Immersed Boundary Methods (IBMs) using weighted least-squares (WLS) [21–23] provide less dispersive and dissipative solutions near the body interface but are computationally intensive and require specific treatment on non-uniform meshes. They also struggle with strong discontinuities, common in hypersonic flows, which can affect the solution and cause crashes.

To overcome these challenges, a low-order interpolation method based on inverse distance weighting is used. This method preserves local maxima, produces smooth reconstructions, and is suitable for industrial configurations due to its applicability on non-uniform meshes and avoidance of negative weighting coefficients. Once variable values at the mirror points are computed, corresponding values at ghost cells are determined using a linear approximation along the boundary normal direction, incorporating prescribed boundary conditions. This

approach ensures appropriate weighting coefficients for neighboring cells, enhancing the stability and accuracy of numerical simulations in the flow solver. Another disadvantage of the ghost-cell-based sharp interface approach is its difficulty in generalizing to arbitrary geometries. For industrial applications, this approach can lack robustness, especially for sharp geometries like cones or slender geometries thinner than the mesh size. In such cases, ghost cells cannot be established for all cells due to insufficient discretization refinement, preventing the presence of a solid cell within them.

Additionally, defining a ghost cell as the cell intersecting the geometry and adjacent to at least one fluid cell leads to issues during the projection of image points. The projection of the image point can be random and dependent on the position of the solid mesh compared to the fluid domain. This randomness can cause solution divergence or poor capture of quantities of interest. In this paper we proposed a new solution involving using “virtual” ghost cells to mitigate these issues, enhancing the robustness and accuracy of the approach for complex industrial geometries.

#### 3.1.1. Interpolation process

For each ghost cell, one correspondent image point $x_i$ is built by projection along the surface normal such as:

$$x_i = x_g + d \cdot n_g \quad (3.1)$$

here d is a fixed distance set to $\sqrt{3}/l$ where l is the size of the solid cell. This distance ensures that the image point never lies within the solid cell, regardless of the configuration.

For each image point $x_i$, the flow field variables are then interpolated from surrounding cells using an interpolation method based on inverse distance weighting [24]. Moreover, the scheme is very efficient in terms of computation time since the coefficients are computed only once and stored during the initialization process.

The interpolation of a field $\Phi$ evaluated at the image point i is defined as:

$$\Phi_i = \frac{\sum_{j=1}^{N} \omega_j(x_i) \Phi_j}{\sum_{j=1}^{N} \omega_j(x_i)} \quad (3.2)$$

where

$$\omega_j(x_i) = \frac{1}{r(x_j, x_i)^p} \quad (3.3)$$

is the weighting function. $x_i$ and $x_j$ are the position of the image point i and the position of the center of the surrounding cell center j, r is the Euclidean distance between these two points, and p a positive integer number usually set to 3 for a 3D resolution. Fig. 3 illustrates the methodology employed. For each image point generated, the interpolation support is built with these neighboring cells contained inside a circle (or sphere in 3D) of fixed size $\sqrt{3}/l$ where l is the size of the cell enclosing the image point. For each node contained inside the circle that is not tagged as a solid node, the corresponding cell is added to interpolation support.

#### 3.1.2. Boundary conditions

Once the flow variable is interpolated at the image point, the boundary condition $\phi_b$ is imposed by applying the interpolated value to the corresponding ghost cell. For Dirichlet condition, the ghost cell value $\phi_g$ is imposed as:

$$\phi_g = 2\phi_b - \phi_i \quad (3.4)$$

In our case, we impose a slip condition on the velocity:

$$\mathbf{u}_g = 2\mathbf{u}_i^t - \mathbf{u}_i \quad (3.5)$$

where $\mathbf{u}_i^t = \mathbf{u}_i - \mathbf{u}_i \cdot \mathbf{n}$ is the tangential velocity.

For Neumann conditions the value of $\phi_g$ is given by

$$\phi_g = \phi_i - \Delta l \left( \frac{\partial \phi}{\partial \mathbf{n}} \right) \quad (3.6)$$
<figure><img src="image_7.png" /><figcaption>Fig. 3. Illustration of the generation of the image points, white circles are the image points, green triangles are the surrounding cells used to create the interpolation support for the image points illustrated by the green circle, red circles are the ghost cells.</figcaption></figure>

where $Δl = |x_i - x_g|$ and $\frac{\partial \phi}{\partial n}$ is the value of gradient at the interface.

In our case, we enforce an adiabatic wall condition (i.e. $\frac{\partial T}{\partial n} = 0$) and zero pressure gradient (i.e. $\frac{\partial p}{\partial n} = 0$). The flow variables imposed on the ghost cells are then:

$$ \left\{ \begin{array}{l} \mathbf{u}_g = 2\mathbf{u}_i^t - \mathbf{u}_i \\ T_g = T_i \\ p_g = p_i \end{array} \right. \qquad (3.7) $$

Once the flow variables are imposed in the ghost cell, the other primitive variables $p_g$ and $E_g$ are computed with (2.3).

### 3.1.3. Improvement for slender geometry

As previously explained, one of the problems with the sharp interface method is the generalization to complex, sharp, or slender geometries. Fig. 4 shows two examples. In these two cases, the image point could be generated in different ways depending on the direction of the projection. The projection of the image point as well as the normal used to impose the slip conditions are thus random and dependent on the position of the solid mesh compared to the fluid domain involving in some cases a divergence of the solution or a poor capture of the quantities of interest.

Our method involves using “virtual” ghost cells, which are generated for each interface between a fluid cell and a cell that has been identified as intersecting the geometry. To better grasp the method, let us consider the example of a thin flat plate (thinner than the Cartesian mesh resolution) (Fig. 5). In this case, multiple ghost cells are generated at identical positions. As a consequence, instead of creating only 5 ghost cells, a total of 11 cells are produced: 5 for the upper geometry faces, 5 for the lower ones, and one additional cell for the left side. These cells are virtual entities and do not have a physical presence within the mesh.

For this purpose, we propose a new approach using “virtual” ghost cells, illustrated in Fig. 5. For each face between a solid and a fluid cell, a “virtual” ghost cell is generated. We detect the nearest surface element by using multiple ray-casts in the direction of the face normal. The normal is then retrieved from the surface element and stored in the data structure for each virtual ghost cell generated.

The challenge now involves determining the specific direction for projecting the image point, which varies for each virtual ghost cell positioned in the same locations. To address this, we use a ray-tracing algorithm, assuming the surface normals of the object are well-defined. For each face of the virtual cell, we define the face using an array of NxN points, as illustrated in Fig. 6. In this scenario, a ghost cell is

intersected by a segment of a surface element (depicted in blue). One of the ghost cell's faces is divided into sub-faces. For every node of these sub-faces, a ray is cast in alignment with the face's normal direction. Some rays intersect with the geometry (shown in red), while others do not (indicated in green). The surface element hit by the ray closest to the face under consideration is subsequently used to determine the direction for projecting the image point. This iterative procedure is applied to all virtual cells defined at identical positions. Fig. 7 provides an illustration of this process in the context of a thin plate and Fig. 8 shows the three image points generated. As we can see, for the ghost cell at the border of the geometry, these points are generated in the direction of the surface normal.

## 3.2. Method 2: Diffuse interface

Let $\phi$ be the solution of the system $\mathcal{N}(\phi)$, where the operator $\mathcal{N}$ defines a generic partial differential equations system. One defines the following notations:

* $\Omega_f$, the fluid domain (Eulerian space)

* $\Omega_S$, the structure domain (Lagrangian space)

* $\Gamma$, the surface of $\Omega_S$, thus the boundary between $\Omega_f$ and $\Omega_S$

* $\mathbf{x} = (x, y, z)$, coordinates in $\Omega_f$

* $\mathbf{q} = (q, r, s)$, coordinates in $\Omega_S$

* $\mathbf{X}(\mathbf{q}, t) = (X(\mathbf{q}, t), Y(\mathbf{q}, t), Z(\mathbf{q}, t))$, the Lagrangian markers describing the surface $\Gamma$

In the rest of the formulation definition, the following convention is adopted: small case letters describe Eulerian variables and capital letters describe Lagrangian variables (here on the boundary $\Gamma$)

The solution of the immersed boundary problem is:

$$ \mathcal{N}(\phi(\mathbf{x}, t)) = f(\phi(\mathbf{x}, t)) \qquad (3.8) $$

$$ \Phi(\mathbf{X}(\mathbf{q}, t)) = \Phi_q = \int_{\Omega_f} \phi(\mathbf{x}, t) \delta(\mathbf{x} - \mathbf{X}(\mathbf{q}, t)) d\mathbf{x} \qquad (3.9) $$

$$ \phi(\mathbf{x}, t) = \int_{\Gamma} \Phi_q \delta(\mathbf{x} - \mathbf{X}(\mathbf{q}, t)) d\mathbf{q} \qquad (3.10) $$

where $f(\phi(\mathbf{x}, t))$ is the volumic IB forcing term imposed at fluid cell close to the solid interface. Eqs. (3.9) and (3.10) are respectively the transfer of the variable $\phi$ from Eulerian space to Lagrangian's, called interpolation step and the transfer from Lagrangian to Eulerian, called spreading step. The information transfer is insured by a convolution with a $\delta$ Dirac function also called interpolation kernel.

In order to simplify the notation, in the following $\mathbf{X}_k = \mathbf{X}(\mathbf{q}_k, t)$ will describe the position of the Lagrangian marker at index $k$ and $\Phi_k$ denote a generic variable defined in the Lagrangian space at the marker $k$. Similarly in the Eulerian space $\Phi_j = \phi(\mathbf{x}_j, t)$ define a general variable at the index $j$.

### 3.2.1. Discretized Dirac function

In the previous part, we presented the continuous form of variable transfer between the two spaces. These transfer methods introduce a Dirac interpolation kernel $\delta$ function and have to be estimated in the discretized space. In the following, we describe the discretized form of the Dirac function for a uniform and a non-uniform Eulerian cell distribution.

**Uniform mesh formulation:** The three-dimensional approximation of the delta function for uniform spatial distribution is given by the product:

$$ \delta_h(\mathbf{x}_j - \mathbf{X}_k) = \frac{1}{H_{X_k}} \delta_h(\frac{x_j - X_k}{H_{X_k}}) \frac{1}{H_{Y_k}} \delta_h(\frac{y_j - Y_k}{H_{Y_k}}) \frac{1}{H_{Z_k}} \delta_h(\frac{z_j - Z_k}{H_{Z_k}}) \quad (3.11) $$

$H_{X_k}$, $H_{Y_k}$, $H_{Z_k}$ are the characteristic cell size in the Lagrangian support $k$ and in a uniform support mesh we have $H_X = H_Y = H_Z = dx$
<figure><img src="image_8.png" /><figcaption>Fig. 4. Illustration of the problem encountered for complex geometries, red circle: ghost point which requires additional processing to generate the image point.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 5. Generation of the virtual ghost cell on a thin flat plate.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 6. Illustration of the ray-tracing process used to detect the direction of the projection of the image point.</figcaption></figure>

where $dx$ is the mesh discretization. The interpolation kernel $\delta_h$ is the discretized form of the Dirac function. In the present work, we use the approximation proposed by Berger and al. [2] :

$$ \delta_h(r) = \begin{cases} \frac{1}{3}(1 + \sqrt{1 - 3r^2}) & 0 \le r \le 0.5 \\ \frac{1}{6}(5 - 3r - \sqrt{1 - 3(1-r)^2}) & 0.5 < r \le 1.5 \\ 0 & \text{otherwise} \end{cases} \quad (3.12) $$

Non-uniform mesh formulation : In 2011, A. Pinelli [4] proposed a generalized window function $\tilde{\delta}$ for non-uniform mesh where the computation procedure to determine $H_{X_k}$, $H_{Y_k}$ and $H_{Z_k}$ is described.

The discrete Dirac function $\delta_h$ is modified using a 2^nd^-order Reproducing Kernel Particle Method. In 3D, the new discrete Dirac function $\tilde{\delta}_h$ is computed as :

$$ \tilde{\delta}_h(\mathbf{x}_j - \mathbf{X}_k) = [b_0 + b_1(\mathbf{x}_j - \mathbf{X}_k) + b_2(\mathbf{y}_j - \mathbf{Y}_k) + b_3(\mathbf{z}_j - \mathbf{Z}_k) $$

$$ \begin{aligned} & + b_4(\mathbf{x}_j - \mathbf{X}_k)(\mathbf{y}_j - \mathbf{Y}_k) + b_5(\mathbf{y}_j - \mathbf{Y}_k)(\mathbf{z}_j - \mathbf{Z}_k) + b_6(\mathbf{z}_j - \mathbf{Z}_k)(\mathbf{x}_j - \mathbf{X}_k) \\ & + b_7(\mathbf{x}_j - \mathbf{X}_k)^2 + b_8(\mathbf{y}_j - \mathbf{Y}_k)^2 + b_9(\mathbf{z}_j - \mathbf{Z}_k)^2] \delta_h(\mathbf{x}_j - \mathbf{X}_k) \end{aligned} $$

Coefficients of the polynomial function $\mathbf{b} = (b_0, ..., b_9)^T$ involve the resolution the linear system :

$$ [m_{a,b,c}] \mathbf{b} = \mathbf{e}_0 \quad (3.13) $$

where

$$ \mathbf{e}_0 = \begin{pmatrix} 1 \\ 0 \\ \vdots \\ 0 \end{pmatrix} \quad (3.14) $$

and

$$ m_{a,b,c} = \sum_{j \in D_k} (x_j - X_k)^a (y_j - Y_k)^b (z_j - Z_k)^c \delta_h(\mathbf{x}_j - \mathbf{X}_k) v_j \quad (3.15) $$

$v_j$ defines the volume of the jth fluid cell, $[m_{a,b,c}]$ is the discrete form of the moment matrix.

### 3.2.2. Interpolation step

The discretized form of the interpolation of a generic variable $\phi$ at the Lagrangian marker $k$ is defined as (see [25]) :

$$ I[\phi]_k = \int_{D_k} \phi(\mathbf{x}, t) \delta(\mathbf{x} - \mathbf{X}(\mathbf{q}, t)_k) d\mathbf{x} \approx \sum_{j \in D_k} \phi_j \delta_h(\mathbf{x}_j - \mathbf{X}_k) v_j \quad (3.16) $$

where $D_k$ is Eulerian support of the Lagrangian marker $k$, $\phi_j$ is the generic fluid quantity at cell j and $v_j$ its volume. $\delta_h$ is the discrete form of the Dirac function either in its uniform mesh formulation $\delta_h$ or non-uniform's $\tilde{\delta}_h$.

### 3.2.3. Determination of Lagrangian quadratures - Spreading step

In this work, two methods for determining Lagrangian marker quadratures are investigated to enforce the boundary condition during the spreading step.

To simplify notations, the following $\mathbf{X}_k = \mathbf{X}(\mathbf{q}_k, t)$ will describe the position of the Lagrangian marker at index k and $\Phi_k$ denote a generic
<figure><img src="image_11.png" /><figcaption>Fig. 7. Generation of the virtual ghost cell on a thin flat plate.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 8. Generation of the images points at the border of the thin flat plate.</figcaption></figure>

variable defined in the Lagrangian space at the marker k. Similarly in the Eulerian space $\phi_j = \phi(x_j, t)$ define a generic variable at the index j.

**Implicit method** : The first method is detailed by J. Favier in [26], where the conservation is assured by the evaluation of a reciprocity parameter. This parameter is obtained by using the conservation property between the Eulerian space and the Lagrangian's. This leads to the resolution of a sparse linear system :

$$A\epsilon = 1 \qquad (3.17)$$

where 1 is the unity vector, $\epsilon$ are the Lagrangian reciprocity parameters. The matrix A is defined by the generic term $A_{kl}$ as follow :

$$A_{kl} = \sum_{j \in D_k} \sum_{l \in D_j} \delta_h(x_j - x_l) \delta_h(x_j - x_k) v_j \qquad (3.18)$$

$A_{kl}$ is defined as the product between the interpolation kernel at the Lagrangian marker k and the interpolation kernel defined in the other Lagrangian markers l. A is diagonal dominant sparse matrix and $A_{kl}$ is not null only when Lagrangian marker l and k are close to each other. The difficulty leads to the determination of the Lagrangian's marker position to prevent an ill-conditioned system and assure conservation during the interpolation-spreading process. As shown in [4,26], the distance between two Lagrangians must be close to the Eulerian cell size and each Eulerian cell that intersects the surface of geometry must contain at most one Lagrangian marker. In this case, the spreading operation is achieved by :

$$S^l(\Phi)_j = \int_{D_j} \Phi_q \delta(x - X(q, t)) dq \approx \sum_{k \in D_j} \Phi_k \delta_h(x_j - X_k) \epsilon_k \qquad (3.19)$$

The exponent notation $S^l$ stands for the implicit determination of the quadrature terms.

**Explicit method** : In 2019, S. Gsell [27] has quantified the reciprocity error $I[S|\Phi]] \approx \kappa \Phi$, where $\kappa = \int \delta(x)^2 dx$, by assuming that the local curvature of the $\Gamma$ surface remains small locally and the variation of $\Phi$ along the boundary are low. This analytical formulation only depends on the Dirac delta function, no linear system has to be solved and, consequently, this reduces the numerical cost compared to usual

implicit methods. In the present work, the discrete value of $\kappa$ at index k reads

$$\kappa_k = \int_{D_k} \delta(x)^2 dx \approx \sum_{j \in D_k} \sum_{l \in D_j} \delta_h(x_j - x_l) \delta_h(x_j - x_k) v_j \qquad (3.20)$$

As for (3.18), the product of Dirac delta function is not null only if markers k and l are close enough to each other. In the article, S. Gsell has shown that one can see the factor $\frac{1}{\kappa}$ as the analytical equivalent of the implicit quadrature term $\epsilon$, previously described. Therefore, the explicit formulation of the spreading is :

$$S^E(\Phi)_j = \int_{D_j} \Phi_q \delta(x - X(q, t)) dq \approx \sum_{k \in D_j} \Phi_k \delta_h(x_j - x_k) \frac{1}{\kappa_k} \qquad (3.21)$$

### 3.2.4. IB forcing term

In the previous section, we presented the general and discrete formulation of the IB method. Here we will describe the general formulation of IB forcing term and how to use it to apply general boundary conditions. When applying (3.8) on a general form of advection equation such as the momentum equation of Euler/Navier-Stokes equations, the boundary conditions are simulated by the forcing term $f(\phi)$ in the Eulerian space and is applied on the right-hand side as follow :

$$\frac{\partial \phi}{\partial t} + \nabla \cdot F(\phi) = f(\phi) \qquad (3.22)$$

where $F(\phi)$ represents the inviscid fluxes in Euler/Navier-Stokes equations. Assuming Eq. (3.22) is satisfied in the Lagrangian space, we can express it explicitly by showing the objective on the field we want to obtain under the considered Lagrangian marker. The procedure employed to evaluate the forcing term in the Lagrangian space $F(\Phi)$ is described in the time-discrete space. In most explicit IBM formulations, the time-stepping procedure uses a predictor-corrector scheme. With a discrete formulation of the temporal term using a first-order scheme, the equation can be expressed as follows:

$$\frac{\Phi_t - \Phi^n}{\Delta t} = -\nabla \cdot F(\Phi^n) + F(\Phi^n) \qquad (3.23)$$

where $\Phi_t$ represents the objective on the field at the Lagrangian maker. The method first performs a predictive step denoted with the suffix *, without considering any solids :

$$\frac{\Phi^* - \Phi^n}{\Delta t} = -\nabla \cdot F(\Phi^n) \qquad (3.24)$$

Substitute Eqs. (3.24) to (3.23) and we get the general formulation of the IB source term in the Lagrangian space :

$$F(\Phi^n) = \frac{\Phi_t - \Phi^*}{\Delta t} \qquad (3.25)$$

Once the IB forcing term is evaluated in the Lagrangian space, its equivalent in Eulerian space at time n is computed with the spreading operator as follows:

$$f(\phi^n) = S[F(\Phi^n)] \qquad (3.26)$$

In the numerical solving of Euler equations, the most classical choices of boundary conditions are slip condition for the velocity i.e the normal
<figure><img src="image_13.png" /><figcaption>Fig. 9. Projection of Lagrangian markers, blue circles represent the cells used to generate wall Lagrangian markers through projection, while black stars indicate the wall Lagrangian markers, and white stars denote the gradient Lagrangian markers. The green box surrounding the green star illustrates the interpolation support used to interpolate fluid quantities at the specific gradient Lagrangian marker.</figcaption></figure>

component of the velocity is set to the solid normal velocity and Neumann condition for the pressure and temperature i.e. zero gradients in the wall-normal direction. The general formulation for these boundary conditions in the Lagrangian space leads to :

$$ \left\{ \begin{array}{l} U \cdot n = 0 \\ \nabla \Phi \cdot n = 0 \end{array} \right. \qquad (3.27) $$

where $\nabla \Phi$ represents the gradient vector of $\Phi$ and $n$ the unit normal.

The novelty of the approach consists of applying the Neumann condition using a second Lagrangian marker projected normally to the surface. Fig. 9 illustrates the process for a set of markers. The generated Lagrangian layer is called the “gradient” layer. In the present case, each Lagrangian marker is projected with a fixed distance $d = \sqrt{3}l$ where $l$ is the size of the corresponding intersecting cell. Using the relations defined in the previous section, the field $\Phi$ is interpolated on the gradient layer and then imposed on the wall markers as follows:

$$ \left\{ \begin{array}{l} \Phi^* = I[(\phi^*)^{(0)}] \\ \Phi_t = I[(\phi^*)^{(1)}] \end{array} \right. \qquad (3.28) $$

where index (0) denote the wall Lagrangian marker and (1) the marker located on the gradient layer.

In the following, we going to present the procedure used to apply the IB forcing term to Euler equations. The boundary conditions to apply lead to :

$$ \left\{ \begin{array}{l} u \cdot n = 0 \\ \nabla T \cdot n = 0 \\ \nabla p \cdot n = 0 \end{array} \right. \qquad (3.29) $$

From (3.22) the IB forcing term is added directly on the right-hand side of Euler equations :

$$ \frac{\partial \mathbf{Q}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{Q}) = f(\mathbf{Q}) \qquad (3.30) $$

The objective is, therefore, to find a formulation of the IB forcing term applied to conservative variables $\rho$, $\rho u$ and $\rho E$ which respects the boundary conditions (3.29).

**IB forcing term on momentum conservation law** : To enforce the slip boundary condition the normal component of the velocity at the Lagrangian point must be equal to zero in the case of a stationary object.

We perform the interpolation of $(\rho u)^*$ on the wall Lagrangian markers:

$$ (\rho U) = I[(\rho u)] \qquad (3.31) $$

From (3.25) and (3.27) applied on the field $\rho u$ at the Lagrangian marker $k$, the IB force can be express as follow :

$$ F((\rho U)^n) = \frac{1}{\Delta t} ((\rho U)^t - (\rho U)^n) = -\frac{1}{\Delta t} ((\rho U_k) \cdot n_k) \qquad (3.32) $$

Then the spreading step is performed to recover the IB force in the Eulerian space. The full expression of $f(\rho u)$ evaluate at time $n$ leads to :

$$ f((\rho U)^n) = S[F((\rho U)^n)] \qquad (3.33) $$

**IB forcing term on energy conservation law** : To simulate the zero gradients on the temperature, we use the perfect gas relation $e_{tot} = e_i + e_u = c_v \cdot T + \frac{1}{2} u^2$ where the kinetic energy is removed to keep the component which only depends on the temperature :

$$ (\rho e_T)^* = (\rho e_{tot})^* - \frac{1}{2} (\rho^*(u^*)^2) \qquad (3.34) $$

The internal energy $(\rho e_i)$ at time * is then interpolated in both the wall and gradient Lagrangian layer where the objective field corresponds to the gradient layer one. From (3.28) :

$$ \left\{ \begin{array}{l} (\rho E_i)^* = I[(\rho e_i)^{(0)}] \\ (\rho E_i)^t = I[(\rho e_i)^{(1)}] \end{array} \right. \qquad (3.35) $$

From (3.25) the IB force leads to :

$$ F((\rho E_i)^*) = \frac{(\rho E_i)^t - (\rho E_i)^*}{\Delta t} \qquad (3.36) $$

From (3.35), (3.36) and including the spreading step, the IB forcing term $f(\rho e)^n$ in the Eulerian space can then be express as :

$$ f(\rho e_i)^n = S[F((\rho E_i)^*)] \qquad (3.37) $$

To take into consideration the term dependent on kinetic energy, the addition of a corrective term is necessary. We propose to use the IB force applied to the momentum equation. The resolution of the momentum equation is resolved before the one on energy and the velocity is updated to recover the velocity at time $n+1$. The corrective term on the energy equation is then calculated as follows:

$$ f(\rho e_u)^n = f(\rho u)^n \cdot u^{n+1} \qquad (3.38) $$

with the index $u$ denote the component of the kinetic energy.

The full expression of the total energy source term $f(\rho e_{tot})$ evaluate at time $n$ leads to :

$$ f(\rho e_{tot})^n = f(\rho e_i)^n + f(\rho e_u)^n \qquad (3.39) $$

**IB forcing term on mass conservation law** : We have just described the procedure to obtain the expression of the IB forcing terms on the momentum and energy equations where the slip condition and the zero temperature gradient are imposed respectively. To simulate the zero pressure gradient we again use the perfect gas relation $p = \rho RT$ and directly impose a zero density gradient. Since the zero gradient on the temperature is set with the IB forcing term $f(\rho e_{tot})$, imposing a zero gradient on the density leads to a zero gradient on the pressure. Similarly to the energy, the density at time * is interpolated on both the wall and gradient Lagrangian layer. From (3.28) :

$$ \left\{ \begin{array}{l} P^* = I[\rho^*]^{(0)} \\ P^t = I[\rho^*]^{(1)} \end{array} \right. \qquad (3.40) $$

With (3.25) the IB force leads to :

$$ F(P^*) = \frac{P^t - P^*}{\Delta t} \qquad (3.41) $$

From (3.40), (3.41) and including the spreading step, the IB forcing term $f(\rho)^n$ in the Eulerian space can then be express as :

$$ f(\rho)^n = S[F(P^*)] \qquad (3.42) $$
**Summary of the method:** The complete method used involves a step-by-step process where the equations are solved in the following order: mass equation, momentum equation, and then energy equation. Between each step, the corresponding primitive variables are updated. The flux is calculated only once during the predictor step and considered constant between each step. From (2.5) the fluxes are computed and the conservative variables are updated at time *:

$$Q_l^* = Q_l^n - \frac{\Delta t}{\gamma_l} \sum_{j=1}^{n=6} \sum_{k=1}^{m} f_{jk}(Q_l^n) \cdot s_{jk} \quad (3.43)$$

The mass equation is then resolved by computing the IB force $f(\rho)$:

$$\rho^{n+1} = \rho^* + \Delta t \cdot S[F(P^*)] \quad (3.44)$$

The update of the $\rho u$ field leads to:

$$(\rho u)^{n+1} = (\rho u)^* - \Delta t \cdot S[F((\rho u)^*)] \quad (3.45)$$

After updating the velocity as: $u_i^{n+1} = \frac{(\rho u)^{n+1}}{\rho^{n+1}}$, the field $\rho e_{tot}$ is updated:

$$(\rho e_{tot})^{n+1} = (\rho e_{tot})^* + \Delta t \cdot S[F((\rho E_i)^*)] + f(\rho u)^n \cdot u^{n+1} \quad (3.46)$$

The complete procedure can be summarized as follows:

1. Predictor step, resolved (2.1) without IB forcing term to get $Q^*$.

2. Resolution of the continuity equation

(a) Interpolate $\rho^*$ on both layer and wall Lagrangian markers (3.40)

(b) Compute $F(P)^n$ (3.41)

(c) Spreading step to get $f(\rho)^n$ (3.42)

(d) Resolve the continuity equation and update the density: $\rho^{n+1}$

3. Resolution of momentum equation

(a) Interpolate $(\rho u)^*$ on wall Lagrangian markers (3.31)

(b) Compute $F((\rho U)_L)^n$ (3.41)

(c) Spread to get $f(\rho u)^n$ (3.42)

(d) Resolve the momentum equation to get $(\rho u)^{n+1}$

(e) Update the velocity $u^{n+1}$ with $(\rho u)^{n+1}$ and $\rho^{n+1}$

4. Resolution of the energy equation

(a) Compute $(\rho e_i)^*$ (3.34)

(b) Interpolate $(\rho e_i)^*$ on both layer and wall Lagrangian markers (3.35)

(c) Compute $F((\rho E_i)^*)$ (3.36)

(d) Spread to get $f(\rho e_i)^n$ (3.37)

(e) Compute the corrective term $f(\rho e_u)^n$ (3.38)

(f) Resolve the energy equation to get $(\rho e_{tot})^{n+1}$ with the full expression of $f(\rho e_{tot})$ (3.39)

(g) Update $e_{tot}^{n+1}$ with $(\rho e_{tot})^{n+1}$ and $\rho^{n+1}$

5. Update the primitive variables $T^{n+1}$ and $p^{n+1}$ with (2.3)

## 4. Results and discussion

As a comparative study, several test cases are proposed. The objective is to evaluate the accuracy of near-wall modeling for Method 1 (Sharp interface) and 2 (Diffuse interface) and how the solver can capture the pressure field since this is the quantity of interest in most industrial applications. Convergence on the integrated forces and on the pressure at the stagnation point will be the focus of these validations. Results are compared with those obtained with high-fidelity body-fitted CFD solvers and some results come from the literature. Moreover, we also compare the obtained results with analytical solutions [28], as typically done in most atmospheric reentry codes to integrate forces

<table><caption>**Table 2**
Free-stream conditions.</caption><tr><td>M</td><td>6.97</td></tr><tr><td>P [Pa]</td><td>126.99</td></tr><tr><td>T [K]</td><td>56.0</td></tr><tr><td>AoA [°]</td><td>0.0</td></tr></table>

and deduce trajectory, see Appendix C for more details. For all the test cases studied, simulations are conducted with a CFL number of 0.8 for both Method 1 and Method 2. Moreover, the CPU times between the two methods are not compared in this study, as the solver used is identical and accounts for the majority of the total CPU time. For both methods, the CPU time required to enforce the boundary conditions represents less than 5% of the total CPU time.

The first test case is a sphere in a hypersonic flow. First, a complete analysis of the proposed Method 2 (Diffuse interface) is performed, where the objective is to determine the best formulation to the specifications of an atmospheric re-entry simulation, i.e. robustness, and precision of the result. This study is followed by a comparison of the two methods with both analytical and high-fidelity body-fitted CFD results. The second test case is an industrial application on a typical reentry capsule of the blunt body type, for which experimental data are available. The final case involves applying the two methods to an academic case of the sharp type, which includes attached shocks: the diamond airfoil.

For the convergence of the solution we define the $L_{inf}$-error norm and the $L_2$-error norm respectively as:

$$L_{inf} = \max_{i=1,N} |\rho_i^{n+1} - \rho_i^n| \quad (4.1)$$

$$L_2 = \sqrt{\sum_{i}^{N} (\rho_i^{n+1} - \rho_i^n)^2} \quad (4.2)$$

where $N$ is the number of total cells in the computational domain.

### 4.1. Hypersonic flow past a sphere

In this part, a verification of parameters used for Method 2 is first presented, where we verify the mesh convergence and observe the difference between the different formulations regarding the conservation property and the penetration error in the case of a hypersonic Euler flow. In a second time, we compare the results obtained with Method 1 and Method 2 with those obtained from Mistral-CFD code [29–31] which already proved his robustness and accuracy on many hypersonic applications.

In the following the $L_\infty$-error norm and the $L_2$-error norm are defined respectively as:

$$\|u - u_{autoref}\|_{\infty} = \max_{i=1,N} \left( \frac{u(i) - u_{autoref}(i)}{u_{autoref}(i)} \right) \quad (4.3)$$

$$\|u - u_{autoref}\|_2 = \sqrt{\frac{\sum_{i}^{N} (u(i) - u_{autoref}(i))^2}{\sum_{i}^{N} u_{autoref}(i)^2}} \quad (4.4)$$

where $N$ is the number of field values.

This study is focused on the 3D sphere case of the diameter of 0.5 m where the initial conditions are imposed as (see Table 2):

Fig. 10 shows the fluid domain for the finest case, where the convergence of the solutions for both methodologies are postponed in Table 3.

#### 4.1.1. Verification on the Diffuse interface method (Method 2)

Fig. 11 presents the mesh convergence of the method, for both the explicit $\kappa$ formulation of the reciprocity parameters and the implicit $\epsilon$ formulation. $\text{reddx}$ represents the mesh size near the wall. This value remains constant for a specific grid level, as the meshes are designed to be uniform near the wall. The convergence study is done by observing
<figure><img src="image_14.png" /><figcaption>Fig. 10. Comparison of computation domain, current work on top, autoreference on bottom.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 11. Convergence on constant spread-interpolated field, $κ$ formulation (left), $ε$ formulation (right), red line: with filtering, green line: without filtering, dotted line: $x$, dashed line: $x^2$.</figcaption></figure>

<table><caption>Table 3
Convergence for the sphere case.</caption><thead><tr><th>Number of cells</th><th>L<sup>method1</sup><sub>2</sub></th><th>L<sup>method1</sup><sub>inf</sub></th><th>L<sup>method2</sup><sub>2</sub></th><th>L<sup>method2</sup><sub>inf</sub></th></tr></thead><tbody><tr><td>295 590</td><td>5.01e<sup>-12</sup></td><td>2.02e<sup>-12</sup></td><td>5.01e<sup>-12</sup></td><td>2.02e<sup>-12</sup></td></tr><tr><td>762 221</td><td>1.17e<sup>-15</sup></td><td>1.22e<sup>-18</sup></td><td>1.17e<sup>-15</sup></td><td>1.22e<sup>-17</sup></td></tr><tr><td>1 508 892</td><td>1.98e<sup>-15</sup></td><td>1.54e<sup>-18</sup></td><td>1.98e<sup>-15</sup></td><td>1.54e<sup>-17</sup></td></tr><tr><td>6 793 989</td><td>8.62e<sup>-16</sup></td><td>2.42e<sup>-18</sup></td><td>8.62e<sup>-16</sup></td><td>2.42e<sup>-18</sup></td></tr></tbody></table>

the $L_2$-norm of the error between the spread-interpolated constant unit scalar field on Lagrangian markers and the constant unit scalar field on the Lagrangian markers, for three mesh sizes. Since the explicit method does not need any pre-processing on the Lagrangian marker, we add the convergence curve in the case of a filtered set of markers. The filtering method is the same as the one used for the implicit formulation (see Eq. (B.1)). One can note that filtering the set of Lagrangian markers does not affect the precision of the reciprocal operation. We also see that the implicit formulation is more precise than the explicit method by 5 orders of magnitude.

In the following, we focus our interest on the respect of the boundary conditions i.e. $u.n = 0$, $\nabla.T = 0$, $\nabla.p = 0$ on the Lagrangian marker. In the following, convergences on the imposed boundary conditions are observed. These observations are made on the converged steady Eulerian flow resolved with the present IBM and time-step notations used in the following are the same as in Section 3.2.4.

For the comparison, once convergence off the steady flow is reached, the interpolation of each field after applying the IB force is done and we compare the convergence on the error between the target value imposed $\Phi^T$ with the value obtained at the end of the process $\Phi^{n+1}$. The comparisons are done with both $ε$ and $κ$ formulations.

The mesh convergence is verified on the projected Lagrangian marker velocity onto the surface normal. Fig. 12(a) shows the $L_2$-norm

of the error between projected velocity and the expected target values equal to zero in case of slip boundary conditions without moving body. For both $κ$ and $ε$ formulations, the convergence is almost of the 2nd order.

In figure 12(b) $L_2$-norm of the error between temperature gradient at time-step $n$ and time-step $n+1$ for both $ε$ and $κ$ formulation is plotted. We found that the convergence is quite similar between $κ$ and $ε$ and followed an order between 1 and 2.

As described in 3.2.4, the boundary condition on the mass equation is imposed by a zero gradient on the density and we pointed out that the zero-gradient density imposed at the wall leads to a zero-gradient pressure in the case of a perfect gas. Figs. 12(c) and 12(d) shows the $L_2$-norm of the error between density gradient and pressure gradient between time-step $n$ and $n+1$. We follow the same conclusions as the previous ones where the convergence using $κ$ and $ε$ are quite similar. We found that the convergence for both pressure and density gradient is between order 1 and 2 and we conclude that the method used to impose a zero-gradient pressure is consistent.

We have just seen that the respect of the boundary conditions on the Lagrangian marker follows the mesh convergence. To conclude this study, we propose to look at the associated mass loss for these two configurations. This study makes it possible to quantify the penetration of the fluid through the geometry. The local mass loss along a surface element is defined as the dot product between $ρu$ and the local normal as :

$$d\dot{m} = \rho u.n.dS \qquad (4.5)$$

For the scaling, we proposed to use the far-field values as :

$$d\dot{m}_0 = \rho_0 u_0.n.dS \qquad (4.6)$$

The integration can then be expressed as follows :

$$\frac{\dot{m}}{m_0} = \int \left( \frac{\rho u}{\rho_0 u_0} \cdot \mathbf{n} \right) . dS \qquad (4.7)$$
<figure><img src="image_16.png" /><figcaption>Fig. 12. Mesh convergence L₂-norm error, κ formulation (blue), ε formulation (green), dashed line: x, dotted line: x².</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 13. Convergence on mass loss on the sphere case for ε and κ formulation.</figcaption></figure>

The discrete formulation is then expressed as :

$$ \frac{\dot{m}}{\dot{m}_0} = \sum_{i=0}^{N} \left( \frac{(\rho \mathbf{u})_i}{\rho_0 \mathbf{u}_0} \cdot \mathbf{n}_i \right) \cdot ds_i \quad (4.8) $$

where suffix i denotes the index of a surface element of the geometry where N is the total number of surface elements.

So we want to look at the conservation of $\frac{\dot{m}}{\dot{m}_0}$ through the surface. The quantity is therefore a volume quantity that we must evaluate at the surface of the geometry. We propose to use a dedicated interpolator where the fluid quantities are interpolated from the center of the fluid cells to the center of the surface triangle of the geometry. For this, we used the same interpolation method defined in 3.1 where $x_j$ and $x_i$ are the position of the center of fluid cell j and the position of the center of the surface element i. The integration is then carried by (4.8).

Fig. 13 shows the convergence of the integrated value for both ε and κ formulations. We note a convergence in mesh for the two formulations. The explicit formulation seems to be more precise. This can be explained by the implication of the filtering method to prevent a negative reciprocity term from reducing the number of points.

<table><caption>Table 4
Convergence for the sphere case.</caption><thead><tr><th>Number of cells</th><th>C<sup>method1</sup><sub>d</sub></th><th>C<sup>method2</sup><sub>d</sub></th><th>C<sup>method1</sup><sub>P<sub>max</sub></sub></th><th>C<sup>method2</sup><sub>P<sub>max</sub></sub></th></tr></thead><tbody><tr><td>295 590</td><td>0.863</td><td>0.892</td><td>1.741</td><td>1.753</td></tr><tr><td>762 221</td><td>0.874</td><td>0.901</td><td>1.769</td><td>1.771</td></tr><tr><td>1 508 892</td><td>0.878</td><td>0.912</td><td>1.802</td><td>1.791</td></tr><tr><td>6 793 989</td><td>0.885</td><td>0.915</td><td>1.823</td><td>1.814</td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>Fig. 14. Convergence on $C_d$ for the sphere case, red line Method 1, green line Method 2, dashed line: x, dotted line: x².</figcaption></figure>

### 4.1.2. Comparison with the high fidelity code

Below we compare the results obtained with those from the Mistral-CFD code. The same space and time discretization schemes are used with similar boundary conditions: slip and adiabatic. Fig. 10 shows a comparison of the domains used, with Mistral-CFD at the bottom and the finest mesh for the current work at the top. Table 4 reports the $C_d$ and $C_p^{max}$ values for the two methodologies, where the mesh convergence is observed. Fig. 14 shows the mesh convergence on $C_d$, where a similar order of convergence is observed for the two methods. A stronger convergence for Method 1 (Sharp interface) is noted for finer meshes.

Plots on the pressure and Mach field are shown in Fig. 15 and compared with high-fidelity code. We can note a major difference
<figure><img src="image_19.png" /><figcaption>Fig. 15. Contour of pressure, current work on top, autoreference on bottom [29], (1) Method 1, (2) Method 2.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 16. Comparison on the shock standoff distance for the finest computation, where $\delta/D$ is the shock standoff distance scaled by the diameter of the sphere, red dot: Method 1, green dot: Method 2, blue dot: autoreference [29].</figcaption></figure>

between the two methodologies. It can be seen that the shock standoff distance is better estimated with Method 1. In Fig. 16 are plotted the shock stand-off distance against the $\theta$ angle. As we can see from Method 1 the shock standoff distance is closer to the autoreference. Moreover, the angle of the bow shock seems to be respected for both methodologies but the shock is shifted for Method 2. These results certainly come from the spreading operation which does not allow to precisely impose the boundary condition at the interface of the object.

In Fig. 17, a contour plot of the $C_p$ on the surface is shown for the finest case and in Fig. 18 the $C_p$ distribution is plotted over the $\theta$ angle. Although the fields in the wake and the shock standoff distance are not respected for Method 2, we can see that both methodologies' pressure coefficient at the walls is well captured. In [32], H.G Hornung et al. performed some simulations on blunted body and sphere cases. They deducted, for the sphere case, that the drag coefficient remained constant in hypersonic configuration by increasing the Mach number. In Table 5 we compare the drag coefficient with that obtained with Mistral-CFD, with analytical solution and result coming from [32]. The drag coefficient is slightly overestimated with Method 2 compared with the body-fitted method. The $C_p^{max}$ is better captured for both methods compared with the analytical solution.

## 4.2. Application on reentry spacecraft: the atmospheric re-entry demonstrator

In this part, we focus our interest on the application of a hypersonic flight vehicle. ESA's Atmospheric Re-entry Demonstrator (ARD) is a

<table><caption>**Table 5**
Comparison of $C_d$ and $C_p^{max}$ on the sphere case.</caption><thead><tr><th></th><th>C<sub>d</sub></th><th>C<sub>p</sub><sup>max</sup></th></tr></thead><tbody><tr><td>Method 1</td><td>0.885</td><td>1.817</td></tr><tr><td>Method 2</td><td>0.915</td><td>1.819</td></tr><tr><td>Mistral-Solver [29]</td><td>0.889</td><td>1.827</td></tr><tr><td>Analytical solution [28]</td><td>0.912</td><td>1.804</td></tr><tr><td>H. G. Hornung [32]</td><td>0.904</td><td>-</td></tr></tbody></table>

<table><caption>**Table 6**
Free-stream conditions.</caption><tr><td>M</td><td>9.92</td></tr><tr><td>P [Pa]</td><td>211.3</td></tr><tr><td>T [K]</td><td>55.7</td></tr><tr><td>AoA [°]</td><td>20.0</td></tr></table>

<table><caption>**Table 7**
Convergence values from coarse to finest with Method 1.</caption><thead><tr><th>Number of cells</th><th>L<sub>2</sub></th><th>L<sub>inf</sub></th><th>C<sub>d</sub></th></tr></thead><tbody><tr><td>311 827</td><td>3.40e<sup>-7</sup></td><td>2.11e<sup>-7</sup></td><td>1.143</td></tr><tr><td>759 112</td><td>4.08e<sup>-8</sup></td><td>5.03e<sup>-8</sup></td><td>1.213</td></tr><tr><td>1 359 461</td><td>5.82e<sup>-9</sup></td><td>7.34e<sup>-9</sup></td><td>1.242</td></tr></tbody></table>

<table><caption>**Table 8**
Convergence values from coarse to finest with Method 2.</caption><thead><tr><th>Number of cells</th><th>L<sub>2</sub></th><th>L<sub>inf</sub></th><th>C<sub>d</sub></th></tr></thead><tbody><tr><td>311 827</td><td>1.30e<sup>-7</sup></td><td>1.10e<sup>-7</sup></td><td>1.151</td></tr><tr><td>759 112</td><td>6.02e<sup>-8</sup></td><td>4.33e<sup>-8</sup></td><td>1.208</td></tr><tr><td>1 359 461</td><td>9.85e<sup>-9</sup></td><td>9.04e<sup>-9</sup></td><td>1.236</td></tr></tbody></table>

3-axis stabilized suborbital reentry vehicle launched from the top of Ariane-503 in 1998. This mission involved the launching of the vehicle into space and its safe recovery in the Pacific Ocean. The ARD has an external diameter of 2.80 m and a height of 2.04 m. The pre-flight experience was carried out by Masurel et al. [33] in the S4 wind tunnel of ONERA. We compare the results with those obtained by the CFD code Lore [34]. The study uses the following far-field conditions with an AoA of 20° (see Table 6).

In this application the solver is coupled with a mesh adaptation method which significantly reduces the number of cells used, better responding to industrial expectations View on the computational domain for the finest case is shown in Fig. 19. For the study, we perform three levels of mesh size, where the convergence of the solution is achieved for both methodologies. Results are shown in Tables 7 and 8.

In Fig. 20 the pressure field in the fluid domain is exhibited. As seen in the previous test case, fields are comparable between both methodologies. $C_p$ contour is shown in Fig. 21 where the stagnation point is
<figure><img src="image_21.png" /><figcaption>Fig. 17. Contour of $C_p$ on the surface, current work on top, autoreference on the bottom, (1) Method 1, (2) Method 2.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 18. Comparison of $C_p$ on the sphere case, red line: Method 1, green line: Method 2, blue line: autoreference [29], dashed line: analytical [28].</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 19. ARD computation domain.</figcaption></figure>

correctly resolved for both methodologies. The integrated values are then compared in Table 9. We see that both $C_d$ and $C_p^{max}$ are close to those obtained from the Lore code and with the analytical ones. Finally, the $C_p$ distribution along the symmetry plane is shown in Fig. 22 and compared at some points with the experimental data. We can observe

<table><caption>Table 9
Comparison of $C_d$ and $C_p^{max}$ for the ARD (finest case).</caption><thead><tr><th></th><th>C<sub>d</sub></th><th>C<sub>p</sub><sup>max</sup></th></tr></thead><tbody><tr><td>Method 1</td><td>1.242</td><td>1.828</td></tr><tr><td>Method 2</td><td>1.236</td><td>1.822</td></tr><tr><td>Lore [34]</td><td>1.227</td><td>1.831</td></tr><tr><td>Analytical [28]</td><td>1.264</td><td>1.827</td></tr></tbody></table>

<table><caption>Table 10
Diamond airfoil: Free-stream conditions.</caption><tr><td>M</td><td>2.0</td></tr><tr><td>P [Pa]</td><td>211.3</td></tr><tr><td>T [K]</td><td>55.7</td></tr><tr><td>AoA [°]</td><td>3.0</td></tr></table>

an overall good agreement on the distribution where the result obtained with Newton's theory seems to overestimate the pressure value on some parts of the vehicle. As we can see the maximum error of the distribution for Method 1 is located close to the edge on the back of the geometry.

### 4.3. Supersonic flow past a diamond airfoil

Here, the academic test case of a diamond airfoil is used to test and analyze the two IBMs in a configuration involving attached shock under supersonic/hypersonic conditions. Fig. 23 describes the geometry. We set the ratio $t_c = t_d/c_d = 0.15$ and the geometry angle $\alpha = 10^\circ$. Free-stream flow characteristics are detailed in Table 10.

To analyze the case, the results of both IBMs are compared to the inviscid shock expansion theory [35,36]. From the front of the diamond, oblique shocks are formed followed by expansion waves at the rear edges of the geometry. The inviscid shock expansion theory allows us to compare the oblique shock angle $\beta$ and flow quantities in different regions (front windward, front leeward, rear windward, and rear leeward). In Fig. 24, regions and oblique shock $\beta$ definition are presented. Front windward and leeward regions are denoted $1d$ and $1u$ respectively and rear windward and leeward regions are denoted as $2d$ and $2u$ respectively. Details on the theory can be found in Appendix D.

Mach distribution obtained for both methodologies are shown in Fig. 25, the sharp interface on the left and diffuse interface on the right. Black lines correspond to the theoretical oblique shock angle. As can be seen, the diffuse interface does not yield satisfactory results. The singularity leads to solution divergence during the spreading operations at both the front and rear corners of the geometry, as depicted in Fig. 25 (right). This divergence is a consequence of the interpolation/spreading process, which is valid only if the interpolated field is sufficiently regular. It highlights the necessity for additional improvement of the diffuse interface when dealing with attached shocks which have not been
<figure><img src="image_24.png" /><figcaption>Fig. 20. Contour of pressure around the ARD surface, (1) Method 1, (2) Method 2.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 21. Contour of $C_p$ on the surface, (1) Method 1, (2) Method 2.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 22. Comparison of pressure field on the ARD, red line: Method 1, green line: Method 2, blue line: autoreference (CFD) [34], dashed line: analytical [28], dark cross: autoreference (experiment) [33].</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 23. Diamond airfoil description.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Fig. 24. Description of observed regions and oblique shock angle $\beta$, diamond airfoil.</figcaption></figure>

carried out in the present work. Conversely, using a sharp interface, as shown in Fig. 25 (left), provides a more effective solution to address this issue. Only results obtained with the sharp interface method can be compared efficiently with the analytical ones.

In Table 11, we compare the oblique shock angle $\beta_{1u}$ and $\beta_{1d}$ with thus obtained analytically (see Appendix D). A good assessment of both windward and leeward oblique shock angles is obtained with the sharp formulation with errors $\le 1\%$. In Table 12 Mach numbers in the four regions are compared to analytical solutions. Errors $\le 2\%$ are observed for all regions, confirming the good predictive capabilities of the sharp formulation.

## 5. Conclusion

In this paper, we propose two new approaches of immersed boundary methods for the resolution of the Euler equations in atmospheric
<figure><img src="image_29.png" /><figcaption>Fig. 25. Contour of Mach number, (1): Method 1, (2): Method 2.</figcaption></figure>

<table><caption>**Table 11**
Comparison on the oblique shock angles $\beta$ with analytical results.</caption><thead><tr><th></th><th>Theory</th><th>Current (Sharp)</th><th>Error (%)</th></tr></thead><tbody><tr><th>β<sub>1d</sub></th><td>36.36</td><td>36.04</td><td>0.880</td></tr><tr><th>β<sub>1u</sub></th><td>42.83</td><td>43.22</td><td>0.911</td></tr></tbody></table>

<table><caption>**Table 12**
Comparison on the Mach-number distribution around the diamond airfoil with analytical results for sharp interface method.</caption><thead><tr><th></th><th>M<sub>1d</sub></th><th>M<sub>1u</sub></th><th>M<sub>2d</sub></th><th>M<sub>2u</sub></th></tr></thead><tbody><tr><th>Theory</th><td>1.522</td><td>1.736</td><td>2.232</td><td>2.489</td></tr><tr><th>Current (sharp)</th><td>1.513</td><td>1.761</td><td>1.836</td><td>2.471</td></tr><tr><th>Error (%)</th><td>0.591</td><td>1.440</td><td>1.837</td><td>0.723</td></tr></tbody></table>

re-entry configuration (hypersonic flows). The first one is based on the sharp interface method and the other on the diffuse interface method. For the sharp interface method, we opt for a ghost cell approach. Low-order interpolation methods using an inverse distance weighting algorithm are adopted to manage stability issues encountered when applied to hypersonic regimes. A methodology using virtual cells is developed to handle complex geometries. For the diffuse interface method, new formulations adapted to compressible Euler equations are developed to better suit the flow regime encountered during atmospheric reentries. In particular, specific formulations of the forcing terms for the mass and energy conservation laws are obtained by adapting the method to handle Neumann conditions. The forcing term in the momentum conservation equation is also modified to impose a slip condition at the interface. As verification, we first look at the impact of the various parameters used in the context of the diffuse interface on a spherical case where we can compare the results obtained for two terms of reciprocity $\kappa$ and $\epsilon$, in terms of convergence on the boundary conditions and fluid penetration. We report significant differences between these two formulations when looking at the interpolation-spreading errors. The implicit formulation is found to be more accurate. However, the effect of the correct imposition of the boundary conditions and its convergence is not significant. A study on the mass loss is performed and leads to the conclusion that the fluid penetration is much lower with the explicit formulation than the implicit one. This can be explained by the fact that for the implicit formulation, the method involves a filtering step where the number of Lagrangian markers is reduced compared with the explicit formulation. In practice, the filtering condition presented to reduce the density of the Lagrangian markers is difficult to generalize to any type of geometry. In a comparative study of the sharp and diffuse interface methods, simulations of hypersonic flow around a sphere and an atmospheric reentry capsule are conducted. We show that both approaches produce results that closely match those obtained using body-fitted methods near the wall. However, the diffuse method provides less resolution in the fields, especially in capturing the shock standoff distance. The final case involves applying both methods to a complex case that includes attached shocks: the diamond airfoil. In this scenario, the diffuse method do not yield converged results. Further developments

to adapt the method to this type of shape are needed. However, by comparing the results obtained with the sharp method to analytical results, we demonstrate that the sharp interface method accurately reproduces the physical behavior with errors of less than 2% for the various parameters studied (angles and Mach fields). Finally, this work has the ambition to adapt the immersed boundary methods for industrial implementation in the context of atmospheric reentries, to handle the high gradients generated by hypersonic flows, manage complex geometries (slender or sharp), and effectively capture strong discontinuities such as attached shocks. For applications involving blunt body geometries, it is found that both IB methods yield closed results in the near-wall region. However, the main novelty lies in applying the IB methods to complex geometries with strong discontinuities. It is also demonstrated that the developed sharp method effectively handles such cases, typical of atmospheric reentry, thus proving its efficiency for industrial applications.

### CRediT authorship contribution statement

**Iko Midani:** Writing – review & editing, Writing – original draft, Validation, Methodology, Investigation. **Valentin Ledermann:** Writing – review & editing, Writing – original draft, Validation, Methodology, Investigation. **Eddy Constant:** Writing – review & editing, Supervision, Methodology, Conceptualization. **Martin Spel:** Writing – review & editing, Supervision, Methodology, Investigation. **Laurent Stainier:** Writing – review & editing, Supervision, Methodology, Investigation. **Pierre Boivin:** Writing – review & editing, Supervision, Methodology, Investigation, Conceptualization. **Julien Favier:** Writing – review & editing, Supervision, Methodology, Conceptualization.

### Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Midani Iko reports financial support was provided by RTech. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

We warmly thank the CNES (Centre National des Etudes Spatiales) and the DGA (Direction Générale de l'Armement) for their financial contributions and interventions that made this study possible.

The authors acknowledge the R. Tech company and the “Association Nationale de la Recherche et de la Technologie” for the CIFRE “Conventions industrielles de formation par la recherche” grant 2020-027 (I. Midani PhD thesis) and for the CIFRE grant 2019-0007 DGA (V. Ledermann PhD thesis). Centre de Calcul Intensif d'Aix-Marseille University is acknowledged for granting access to its high performance computing resources.
<figure><img src="image_30.png" /><figcaption>Fig. A.26. Illustration of intersecting cells identification. On left, raytracing process used to tagged the nodes. White circle is a node tagged as outside where the black one is tagged as inside, the intersection are illustrate with red cross. On right, the black triangle is a cell tagged as inside, white triangle a cell tagged at outside and the red circle indicates a cell tagged as intersecting the geometry.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Fig. A.27. Illustration of fluid cell identification, white triangle: fluid cell marked as outside the solid, dark triangle: fluid cell marked as inside the solid, blue circle: cell marked as intersecting cell.</figcaption></figure>

## Appendix A. IBM data structure

For both methodologies, the data structure is constructed in similar ways. This requires the detection of fluid cells that intersect the surface as well as the differentiation between cells internal and external to the geometry. Thus, being D the fluid domain, this is divided into three portions: $\Omega_{solid}$, $\Omega_{inside}$ and $\Omega_{outside}$ and the information is stored inside a marker variable $\xi$ defined as :

$$ \xi(x) = \begin{cases} 0 & \text{if } x = (x, y, z) \in \Omega_{\text{inside}} \\ 1 & \text{if } x = (x, y, z) \in \Omega_{\text{outside}} \\ 2 & \text{if } x = (x, y, z) \in \Omega_{\text{solid}} \end{cases} \quad (A.1) $$

The tagging procedure is based on ray-tracing algorithm. A ray that originates from the location to be checked is considered and the intersections between this ray and the given surface are computed. Fig. A.26 illustrates the method. The origins of the rays are the mesh nodes. For each node within the octree domain that falls inside the bounding box of the object, N rays are cast in random directions, with N = 5 in our case. The number of intersections between each ray and the geometry is then counted. If the count is even, the node is tagged as outside the geometry otherwise the node is tagged as inside. A cell is marked as an intersecting cell if at least one of its octree nodes is classified as outside and another as inside. If all nodes are outside or all are inside, the cell is marked as outside or inside, respectively. The intersection between a ray and the surface is carried out using the geometrical algorithms reported in [37].

Fig. A.27 illustrates the result.

## Appendix B. Verification of the reciprocity parameters

This section will discuss a study on the Lagrangian reciprocity parameters computations presented in part 3.2.3. As J. Favier pointed out in [26], when this implicit method is used, the Lagrangian distribution must be as uniform as possible and at most one marker per Eulerian cell. In our case, we want a general algorithm that manage all kind of geometry and in particular spacecraft shapes, and for that we cannot set the Lagrangian marker by hand. The original markers must be filtered in order to prevent wrong conditioning of the linear system (3.17). The filtering condition is an empirical formulation and in the present algorithm a Lagrangian marker is set as valid (i.e. not filtered) if it respect the following condition for each valid markers within its neighborhood :

$$ \gamma_{\text{fil}} \frac{\| \vec{L} \cdot \vec{n}_c \|}{\| \vec{n}_c \|} \le \frac{\| \vec{n}_L \cdot \vec{n}_c \|}{\| \vec{n}_c \|} \quad (B.1) $$

where $\vec{n}_L$ is the vector between the compared Lagrangian markers, $\vec{n}_c$ is the vector between the cell centers of the Lagrangian markers owner cells (see figure below), $\vec{L}$ is a vector based on owner cells size and it is computed as:

$$ \vec{L} = \begin{pmatrix} \frac{L_{1x}}{2} + \frac{L_{2x}}{2} \\ \frac{L_{1y}}{2} + \frac{L_{2y}}{2} \\ \frac{L_{1z}}{2} + \frac{L_{2z}}{2} \end{pmatrix} $$

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

The following study will try to determine an appropriate value of $\gamma_{fil}$. When the Lagrangian distribution is denser than the Eulerian distribution some of the implicit reciprocity parameters $\epsilon$ obtained from (3.17) become negative and it can imply divergence of the code, since it creates a negative IB source term. Figs. B.28 and B.29 present the effect of the filtering factor on the number of negative $\epsilon$ for $0.6 \le \gamma_{fil} \le 1.4$ and around $\gamma_{fil} = 1$. One can observe the difference in behavior according to the interpolation kernel used.

Now we have ensured the positiveness of the $\epsilon$ reciprocity parameters by choosing an appropriate $\gamma_{fil}$, an investigation on the convergence of the interpolation-spreading operation against the solver tolerance is presented. The sparse linear system (3.17) is resolved using an iterative conjugate gradient algorithm. Fig. B.30 shows the evolution of the $L_\infty$-error norm against the conjugate gradient solver residual.
<figure><img src="image_33.png" /><figcaption>Fig. B.28. Percentage of negative $\epsilon$ vs. $\gamma_{fl}$ for $\delta_h/\delta_h^i$ interpolation kernel.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Fig. B.29. Percentage of negative $\epsilon$ vs. $\gamma_{fl}$ for $\delta_h/\delta_h^i$ interpolation kernel, zoomed views.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Fig. B.30. $L_\infty$ interpolation-spreading error norm against the solver residue for $\delta_h/\delta_h^i$ interpolation kernel.</figcaption></figure>

## Appendix C. Modified Newton theory

The formulation of the modified Newton pressure coefficient at the i face of the surface $C_p^i$ is defined as:

$$C_p^{\text{max}} = \frac{2}{\gamma M^2} \cdot \left( \left( \left( \frac{(\gamma + 1)^2 M^2}{4\gamma M^2 - 2(\gamma - 1)} \right)^{\frac{\gamma}{\gamma-1}} \cdot \frac{1 - \gamma + 2\gamma M^2}{\gamma + 1} \right) - 1 \right) \quad (C.1)$$

$$C_p^i = \begin{cases} C_p^{\text{max}} \cos \alpha^2 & \cos \alpha > 0 \\ 0 & \text{otherwise} \end{cases} \quad (C.2)$$

## Appendix D. Diamond airfoil: analytical solution

The top and bottom angles of the oblique shocks are given by the following equations:

$$\tan(\alpha - \mathcal{A}) = 2 \cot \beta_{1u} \frac{M_0^2 \sin^2 \beta_{1u} - 1}{M_0^2(\gamma + \cos 2\beta_{1u} + 2)} \quad (D.1)$$

$$\tan(\alpha + \mathcal{A}) = 2 \cot \beta_{1d} \frac{M_0^2 \sin^2 \beta_{1d} - 1}{M_0^2(\gamma + \cos(2\beta_{1d}) + 2)} \quad (D.2)$$

Once the shock angle is found, we can express the Mach number after the shocks as:

$$M_{1u} = \frac{1}{\sin^2(\beta_{1u} - \alpha + \mathcal{A})} \frac{1 + \frac{\gamma+1}{2} M_0^2 \sin^2 \beta_{1u}}{\gamma M_0^2 \sin^2 \beta_{1u} - \frac{\gamma-1}{2}} \quad (D.3)$$

$$M_{1d} = \frac{1}{\sin^2(\beta_{1d} - \alpha - \mathcal{A})} \frac{1 + \frac{\gamma+1}{2} M_0^2 \sin^2 \beta_{1d}}{\gamma M_0^2 \sin^2 \beta_{1d} - \frac{\gamma-1}{2}} \quad (D.4)$$

In order to find the post-expansion distribution, we need to evaluate the Prandtl-Meyer function:

$$v(M_{1u/d}) = \sqrt{\frac{\gamma+1}{\gamma-1}} \arctan \sqrt{\frac{\gamma-1}{\gamma+1}} (M_{1u/d}^2 - 1) - \arctan \sqrt{Ma_{1u/d}} - 1 \quad (D.5)$$

Then the Mach number after the expansion can be evaluated with:

$$v(M_{1u/d}) + 2\alpha = \sqrt{\frac{\gamma+1}{\gamma-1}} \arctan \sqrt{\frac{\gamma-1}{\gamma+1}} (M_{2u/d}^2 - 1) - \arctan \sqrt{Ma_{2u/d}} - 1 \quad (D.6)$$

## Data availability

Data will be made available on request.

## References

[1] Peskin CS. Flow patterns around heart valves: A numerical method. J Comput Phys 1972. URL https://doi.org/10.1016/0021-9991(72)90065-4.

[2] Roma AM, Peskin CS, Berger MJ. An adaptive version of the immersed boundary method. J Comput Phys 1999. URL https://url.org/10.1006/jcph.1999.6293.

[3] Iaccarino G, Kalitzin G, Khalighi B. Towards an immersed boundary RANS flow solver. Aerosp Sci Meet Exhibit 2003. URL.

[4] Pinelli A, Naqavi IZ, Piomelli U, Favier J. Immersed-boundary methods for general finite-difference and finite-volume Navier-Stokes solvers. J Comput Phys 2010. URL.

[5] Constant E, Favier J, Meldi M. An immersed boundary method in openfoam: Verification and validation. Comput Fluids 2017;157:55-72. http://dx.doi.org/10.1016/j.compfluid.2017.08.001.

[6] Riahi H, Meldi M, Favier J, Serre E, Goncalves E. A pressure-corrected immersed boundary method for the numerical simulation of compressible flows. J Comput Phys 2018. URL https://url.org/10.1016/j.jcp.2018.07.033.

[7] Kumar V, Sharma A, Singh RK. Central upwind scheme based immersed boundary method for compressible flows around complex geometries. Comput Fluids 2020. URL https://url.org/10.1016/j.compfluid.2019.104349.

[8] Kim J, Kim D, Choi H. An immersed-boundary finite-volume method for simulations of flow in complex geometries. J Comput Phys 2001. URL https://url.org/10.1006/jcph.2001.6778.

[9] Ren W, Shu C, Yang W. An efficient immersed boundary method for thermal flow problems with heat flux boundary conditions. Int J Heat Mass Transfer 2013. URL https://url.org/10.1016/j.ijheatmasstransfer.2013.05.020.

[10] Park SG, Chang CB, Kim B, Sung HJ. Simulation of fluid flexible body interaction with heat transfer. Int J Heat Mass Transfer 2017. URL https://url.org/10.1016/j.ijheatmasstransfer.2017.03.012.

[11] Schneiders L, Günther C, Meinke M. An efficient conservative cut-cell method for rigid bodies interacting with viscous compressible flows. J Comput Phys 2016;311:62-86. http://dx.doi.org/10.1016/j.jcp.2016.01.026.

[12] 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 2018;357:230-62. http://dx.doi.org/10.1016/j.jcp.2017.12.030.
[13] 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 2017;15(01):1750080. http://dx.doi.org/10.1142/s0219876217500803.

[14] Qu Y, Shi R, Batra R-C. An immersed boundary formulation for simulating high-speed compressible viscous flows with moving solids. J Comput Phys 2018;354:672–91. http://dx.doi.org/10.1016/j.jcp.2017.10.045.

[15] Khalili ME, Larsson M, Müller B. Immersed boundary method for viscous compressible flows around moving bodies. Comput Fluids 2018;170:77–92. http://dx.doi.org/10.1016/j.compfluid.2018.04.033.

[16] Toro Eleuterio F. Riemann solvers and numerical methods for fluid dynamics. 1999.

[17] Sod GA. A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws. J Comput Phys 1978. URL https://doi.org/10.1016/0021-9991(78)90023-2.

[18] Péron S, Benoit C, Renaud T. An immersed boundary method on Cartesian adaptive grids for the simulation of compressible flows around arbitrary geometries. Eng Comput 2020.

[19] Cui Z, Yang Z, Jiang H. A sharp-interface immersed boundary method for simulating incompressible flows with arbitrarily deforming smooth boundaries. Int J Comput Methods 2018.

[20] Tran P, Plourde F. Computing compressible internal flows by means of an immersed boundary method. Comput Fluids 2014.

[21] Choung H, Saravanan V, Lee S. Nonlinear weighting process in ghost-cell immersed boundary methods for compressible flow. J Comput Phys 2021;433:110198.

[22] Khalili ME, Larsson M, Müller B. High-order ghost-point immersed boundary method for viscous compressible flows based on summation-by-parts operators. Internat J Numer Methods Fluids 2018;89(7):256–82.

[23] Mohammadi MH, Sotiropoulos F, Brinkerhoff J. Moving least squares reconstruction for sharp interface immersed boundary methods. Internat J Numer Methods Fluids 2019;90(2):57–80.

[24] Franke R. Scattered data interpolation: tests of some methods. Math Comp 1982. URL https://doi.org/10.2307/2007474.

[25] Li Z, Favier J, D'Ortona U, Poncet S. An immersed boundary-lattice Boltzmann method for single and multi-component fluid flows. J Comput Phys 2016. URL https://doi.org/10.1016/j.jcp.2015.10.026.

[26] Favier J. Modélisation et simulation numérique d'écoulements au voisinage de paroisbiologiques ciliées et de revêtements aéronautiques bioinspirés. 2016.

[27] Gsell S, D'Ortona U, Favier J. Explicit and viscosity-independent immersed-boundary scheme for the lattice Boltzmann method.

[28] Andeerson JD. Hypersonic and high-temperature gas dynamics. 1989, p. 64–5.

[29] Grossir G, Puorto D, Ilich Z, Paris S, Chazot O, Rumeau S, Spel M, An- naloro J. Aerodynamic characterization of space debris in the VKI longshot hypersonic tunnel using a free-flight measurement technique. Int Conf Flight Veh Aerothermodyn Re- Entry Mission Eng 2019.

[30] Scott T, Dieudonné W, Spel M. 34Th aiaa fluid dynamics conference and exhibit. 2012, URL https://doi.org/10.2514/6.2004-2423.

[31] Vérant J, Spel M, Bonnal B. EXOMARS capsule aerodynamics analysis. In: 10th AIAA/ASME joint thermophysics and heat transfer conference. 2010. URL https://doi.org/10.2514/6.2010-4646.

[32] Hornung HG, Schramm J-M, Hannemann K. Hypersonic flow over spherically blunted cone capsules for atmospheric entry. Part 1. The sharp cone and the sphere. J Fluid Mech 2019. URL https://doi.org/10.1017/jfm.2019.342.

[33] J. M, Jobard G, Bouis X. Heat flux study and pressure measurements on ard and soyuz vehicles in perfect gas condition in the s4ma wind tunnel. Technical report PV 4/8015 GY 400, ONERA; 2001.

[34] Muylaert J, Walpot L, Vennemann D. A review of European code-validation studies in high-enthalpy flow. In: 20th AIAA advanced measurement and ground testing technology conference. 1999. URL https://doi.org/10.1098/rsta.1999.0430.

[35] Liepmann HW, Roshko A. Elements of gasdynamics. 1957.

[36] Frapolli N, Chikatamaria SS, Karlin IV. Entropic lattice Boltzmann model for gas dynamics: Theory, boundary conditions, and implementation. Phys Rev E 2016. URL.

[37] Intel embree: a high-performance ray tracing library. 2023. https://www.embree.org/. [Accessed 17 March 2023].