Assessment of Immersed Boundary Methods for Hypersonic Flows with Gas-Surface Interactions

Immersed boundary (IB) methods with adaptive mesh refinement (AMR) techniques are assessed for atmospheric entry applications, including effects of chemical nonequilibrium (CNE) and gas-surface interactions (GSI). The performance of a conservative cut-cell and two non-conservative ghost-cell IB methods is assessed in comparison with analytical solutions, data from literature, and results obtained with a reference solver that operates on body-fitted grids. All solvers use the same external thermochemistry library so that all observed differences can be attributed to the underlying numerical methods. Results from eight benchmark cases are reported. Four cases are selected to verify the implementation of chemistry, transport properties, catalytic boundary conditions, and shock capturing. Four validation cases consider blunt geometries with adiabatic/isothermal and inert/catalytic/ablative boundary conditions. Overall, the results obtained with the IB solvers are in very good agreement with the reference data. Discrepancies arise with ghost-cell methods for cases with large temperature or concentration gradients at the wall and are attributed to mass conservation errors. Only a strictly conservative cut-cell IB method is on par with body-fitted grid methods.


Introduction
Hypersonic flows experienced during atmospheric entry of capsules or space debris are characterized by strong shock waves and thermochemical nonequilibrium effects through the excitation of the internal energy modes of species and rapid chemical reactions in the shock layer.The hot gas interacts with the surface thermal protection system (TPS) material installed to protect the spacecraft from this hostile environment.Depending on the characteristics of the TPS material, these gas-surface interactions (GSI) involve catalysis as well as ablation.While the former accelerates the exothermic recombination reactions leading to increased heat transfer towards the surface, the latter alleviates the heat load by means of physicochemical decomposition and mass loss.These ablative GSI change the shape of the object by surface recession.Understanding these interactions is crucial for predicting the surface stresses and heat fluxes, as well as the uncontrolled trajectory of space debris.Ground testing is indispensable for validation purposes; however, no facility can simultaneously replicate all aspects of atmospheric entry flows [1].Hence, computational fluid dynamics (CFD) simulations are essential for the efficient aerothermodynamic analysis and design of future spacecraft.
Most CFD solvers used for high-speed and high-enthalpy applications employ body-fitted structured grids [2,3,4].In these solvers, alignment of the grid with the shock and the surface needs to be ensured for an accurate prediction of the flow field.Generating these types of grids usually involves strenuous effort from the user especially for detailed features and incremental geometry updates [5].Unstructured grids have also been explored; however, issues affecting the heat flux predictions at the surface were reported [6,7].
A promising alternative is the use of adaptive mesh refinement (AMR) techniques based on piecewise Cartesian grids with immersed boundary (IB) methods.There has been a recently growing interest in IB-AMR solvers for atmospheric entry applications [8,9,10,11,12], mainly for their potential in considering complex and deforming geometries, and better robustness and higher computational efficiency compared to body-fitted mesh-deformation methods.These methods also allow for a relatively straightforward implementation of high-order schemes.However, special care must be taken to have sufficient grid resolution near the boundaries, as it is more difficult for immersed boundary methods to efficiently resolve thin boundary layers over curved surfaces.To address this shortcoming, a blend of Cartesian grids in the fluid and body-fitted grids near the surface can be employed [10,11].This approach has been successful in reducing the required number of cells and providing better resolution of the thermal boundary layer.In general, a blended grid approach is well suited for shapes with smooth curvatures.However, it is susceptible to the same drawbacks inherent to body-fitted grids, for instance, their difficult adaptation to complex deforming geometries.
Arslanbekov et al. [8], Sekhar and Ruffin [9], and more recently Brahmachary et al. [12] demonstrated the benefits of using IB-AMR solvers for a number of relevant cases.These studies have generally indicated good predictions for wall pressure and skin friction distributions, while emphasizing the difficulty in accurately predicting wall heat fluxes.As with more recent contributions [10,11], these studies were mostly performed with ghost-cell methods from the family of discrete forcing IB approaches [13].Ghost-cell methods impose boundary conditions by extrapolating the fluid solution into the solid, i.e. into ghost cells that are neighbouring fluid cells.Relying solely on ghost-point extrapolation does not ensure strict conservation of mass, momentum, and energy.A strictly conservative approach is the cut-cell finite-volume method, which splits fluid and solid domains into consistently deformed finite volumes.The implementation of a cut-cell method for three dimensions and high-order schemes is not as straightforward as the ghost-cell approach, and it also introduces additional challenges such as cut-cells with very small fluid volumes.However, the main advantage of the cut-cell method lies in satisfying the conservation of mass, momentum, and energy near the wall [13].
In this paper, ghost-cell and cut-cell IB methods are scrutinized through a curated list of benchmark case studies relevant for atmospheric entry applications.Main contributions of this paper are twofold: • to assess the accuracy of IB methods for applications with strong thermal gradients and gas-surface interactions.
• to establish a set of well-defined test cases for the verification and validation of IB methods for atmospheric entry.
Results obtained with the IB-AMR solvers INCA [14,15] and CHESS [16,17] are compared to reference results obtained with the body-fitted finite-volume solver US3D [6] in addition to data from literature.A consistent comparison to assess the accuracy of the numerical methods is achieved by coupling each of the flow solvers with the same thermochemistry library, Mutation ++ [18].
The paper is structured as follows: Governing equations and modelling approaches are presented in Section 2. Solver methodologies are introduced in Section 3. Results of the benchmark case studies are presented and discussed in Section 4, while the influence of the different IB methologies is further investigated in Section 5. Concluding remarks are made in Section 6.

Governing Equations
The compressible Navier-Stokes equations are solved in their conservative form for a reacting multicomponent fluid, where ρ i is the species partial density for the i th species, u is the mixture average velocity, ωi is the source term associated with the production or consumption of species due to chemical reactions, ρ is the mixture density, p is the mixture pressure, and E = e + u 2 /2 is the specific total energy, which is the sum of the thermodynamic internal energy e and the kinetic energy.
External forces due to gravitational or electromagnetic effects, and radiative energy exchanges are not considered for the cases in this study.Both solvers considered in this work can perform under thermal nonequilibrium with multi-temperature methods, such as that of Park [19].However, results presented in this paper are obtained with a thermal equilibrium assumption.The ideal gas assumption leads to the equation of state p = ρRT , where R = R/M is the mixture gas constant obtained from the universal gas constant R and the mixture average molar mass M .These mixture properties are modelled according to Dalton's law through their constituent species as Two models are considered for the species diffusion flux J i : Fick's law with a correction to ensure conservation of mass as with the mixture-averaged diffusion coefficients , obtained by Wilke's average of the binary diffusion coefficients D ij .The second diffusion model uses the solution of the Stefan-Maxwell equations, where x i are the mole fractions, and M i are the species molar masses.This formulation is computationally costlier, but theoretically more accurate [20].
Viscosity and thermal conductivity are obtained through a linear system solution using an LDL T decomposition as opposed to the common approach of using simplified mixture rules [21,22,23].The viscous stress tensor τ is defined assuming Stokes' hypothesis as where µ is the dynamic (shear) viscosity of the mixture.The total heat flux vector q includes the contributions from conduction and mass diffusion, where T is the temperature.The first term stems from Fourier's law with the thermal conductivity λ of the mixture, and the second term accounts for the transport of enthalpy by species diffusion, with h i as the species enthalpy.

Physicochemical Modelling
The models used in state-of-the-art CFD solvers capable of simulating the aforementioned phenomena vary considerably.Broadly, choices need to be made on the thermodynamic database, the treatment of TCNE effects, the transport properties modelling, and the approach for handling GSI.For details we refer to several published studies that evaluate the impact of these selections in modelling thermal nonequilibrium [24,25], species diffusion [20,26], viscosity and thermal conductivity [26,27], rate of catalysis [28], and ablation [29].As important quantities of interest, such as surface heat fluxes, are highly sensitive to the modelling approaches, large discrepancies between the results obtained with hypersonic CFD codes are common [3,4].Hence, this variety of approaches often obscures a clear assessment of the underlying numerical methods, when comparing different solvers.
Based on these considerations, each of the solvers used in this study is coupled with the multicomponent thermodynamic and transport properties for ionized gases in C++ (Mutation ++ ) open-source library.Mutation ++ provides all required physiochemical models for thermodynamics, transport properties, chemical kinetics, and GSI.A detailed description of Mutation ++ is presented by Scoggins et al. [18].
The caloric properties of the species are approximated with standard NASA-9 polynomial fits [30].Species mass diffusivities, viscosities, and thermal conductivities are provided by Mutation ++ according to multicomponent Chapman-Enskog formulations [31].The chemical reaction mechanisms, that is, species mass rates and their analytical Jacobians with respect to species densities, are also provided by Mutation ++ .
Catalytic and ablative surface boundary conditions are imposed by solving a mass balance [32,18], with v blow as the surface-normal blowing velocity, which is nonzero only for an ablative boundary.Terms from left to right refer to convective flux due to blowing, diffusive flux, and species source term due to surface reactions.A probability based approach is employed for computing this chemical source term for the surface, written as ωi,wall = γm i F i,impin , where γ = F i,react /F i,impin is the ratio of reacting to impinging species fluxes and it describes the efficiency of the process, and m i is the mass of the i th species [32].Assuming the species at the wall have a Maxwellian distribution function, the impinging species flux is where k B is the Boltzmann constant, T w is the wall temperature, and n i is the number density of the i th species [33].From the mass blowing rate ṁ = i ωi,wall , the blowing speed is calculated by Values obtained for species densities and mass blowing speeds are then imposed as boundary conditions for the Navier-Stokes equations.

Numerical Methods
We consider three different methodologies for imposing surface boundary conditions in the framework of finite-volume methods.Schematics of the body-fitted, ghost-cell IB, and cut-cell IB approaches are shown in Fig. 1.The arc near the middle of each sketch indicates the surface that demarcates the fluid above from the solid below it.The other black lines are grid lines and the filled dots indicate cell centers.The color code matches the one used for presenting the results in Section 4.
The classical body-fitted grid method, Fig. 1a, simply makes use of the grid's alignment with the geometry.An example stencil is drawn in the sketch and the boundary intercept is indicated by a colored hollow circle.Hollow circles indicate the stencil of the discretization scheme.The fluid-cell solutions and boundary conditions are used to reconstruct quantities at cell interfaces according to the chosen numerical scheme.
The two other approaches use IB methods on Cartesian grids.The ghostcell IB approach [34] seen in Fig. 1b imposes boundary conditions by extending the fluid solution to ghost-cells.The virtual flow solution of ghost-cells is set by extrapolating the nearby fluid solution according to boundary conditions at the nearest fluid-solid interface.Since the grid does not conform with the geometry, the solution at an image point in the fluid needs to be found through interpolation using the surrounding fluid-cell solutions.An example stencil is colored in the sketch, with the fluid points used in the interpolation connected by dotted lines.Ghost-cells are marked with a striped pattern in the sketch.While being relatively straightforward to implement, this approach does not ensure strict conservation of mass, momentum, and energy at the interface between the fluid and the solid.Fluxes are reconstructed from the fluid-cell and the ghost-cell solutions on the Cartesian grid without considering the location and shape of the fluid-solid interface.Errors in implicitly satisfying the conservative flux boundary condition therefore result from the nonlinearity of the flux function, from the image point interpolation, and from the interface curvature.
The cut-cell IB approach [35,36,37], see Fig. 1c, ensures strict conservation by considering the flux balance for the part of the cell that belongs to the fluid domain.These consistently deformed finite volumes and their cell faces are colored in the sketch.Fluxes over the cell faces of the cut cells are scaled according to the wetted areas.The exchange of mass (e.g. with surface reactions), momentum, and energy through the fluid-solid interface is calculated from the prescribed boundary conditions and the local fluid solution.The latter is acquired by interpolation from the surrounding cell values and the boundary conditions.An example stencil is colored in the sketch for the cut-cell interpolation.The other stencil in the sketch is identical to the ghost-cell IB approach.This addition to the cut-cell method refers to the specific implementation within the INCA solver and will be discussed in Section 3.2.Cut-cells with a very small fluid volume fraction require a special treatment to ensure stable time integration.They are typically mixed or merged with nearby cells [37].

Body-fitted Solver
The body-fitted solver considered in this study is US3D, which is a highfidelity flow solver specifically designed for aerodynamic applications in the hypersonic regime by the University of Minnesota and NASA [6].It solves the compressible chemically reacting Navier-Stokes equations in a finite-volume framework on unstructured body-fitted grids.Among the several numerical schemes available in the solver, all simulations carried out within this work use the modified Steger-Warming scheme [38], which is suitable for steady computations.A MUSCL approach [39] is employed to obtain second-order accurate fluxes.Both explicit and implicit time integration methods are available; in this work, rapid convergence to steady state is achieved with the data parallel line relaxation (DPLR) method [40].US3D is equipped with chemistry/multi-temperature source terms and transport properties with the possibility to account for high temperature and high pressure effects.Native routines can be further extended by user-defined subroutines, which allow coupling the solver to external libraries; we refer to Capriati et al. [41] for the coupling with Mutation ++ .

Immersed Boundary Solvers
Two IB solvers are considered: one able to use both the cut-cell and the ghost-cell methods, and another using only the latter.
Employing a cut-cell IB methodology, INCA is a high-fidelity finitevolume solver for direct numerical simulations (DNS) and large eddy simulations (LES) of the compressible chemically reacting Navier-Stokes equations and provides a large number of different discretization schemes on three-dimensional block-Cartesian AMR grids [14,15].For the purposes of this study, a third-order weighted essentially non-oscillatory (WENO) scheme [42] with HLLC flux function [43] is selected to discretize the inviscid terms.WENO schemes permit high accuracy in smooth regions, while ensuring stable and sharp capturing of discontinuities.Second-order centered differences are used for the viscous terms and the explicit third-order Runge-Kutta scheme of Gottlieb and Shu [44] is selected for time integration.Chemical source terms are treated using Strang's second-order time splitting scheme [45] to alleviate the numerical stiffness caused by these terms.The chemical source terms thus reduce to a system of ordinary differential equations, which is solved by the VODE library [46].INCA employs a unique improvement to the common cut-cell methodology [37], which we refer to as the cut-element method [47,48].This method represents the fluid-solid interface through cut-elements, which are derived from the Cartesian mesh and the triangulation of the surface geometry.Instead of considering a planar intersection of a finite-volume cell with the wall surface, as typically done in cut-cell methods [49,50], cut-elements maintain all details of the intersection of the grid with the surface triangulation.The interface within each cut-cell is thus represented by several cut-elements belonging to different surface triangles to yield sub-cell accuracy and robustness for complex geometries.This method is a consistent and conservative extension of the finite volume flux balance to accommodate cells being split by boundaries.Further details on this cut-element methodology and its extension to incorporate GSI and the effects of thermal nonequilibrium are provided in Ref. [51].
INCA employs ghost-cells to allow for the use of unmodified stencils throughout the domain, as shown in Fig. 1c.Moreover, the cut-cell procedure that ensures strict conservation can be switched off to use only the ghost-cell method.We will discuss results obtained with the INCA ghost-cell method for selected cases in Section 5.
In contrast to INCA, the IB method implemented in the flow solver CHESS of Politecnico di Bari [17] fully relies on ghost-cells.The numerical method utilized by CHESS is based on the flux vector splitting proposed by Steger and Warming [38] with a second-order MUSCL reconstruction in space [39] for the hyperbolic terms.Discretization of the viscous fluxes uses Gauss's theorem in conjunction with a second-order linear reconstruction of the solution.A third-order explicit Runge-Kutta scheme is employed for time integration of the transport terms in the Navier-Stokes equations.Following the Runge-Kutta time step, chemical source terms are computed by means of a Gauss-Seidel scheme.CHESS also uses AMR to provide appropriate resolution of shocks and boundary layers [52] and uses the same physicochemical models as US3D and INCA [16].Further details on the solver can be found in the aforementioned works [16,17].

Benchmark Cases
We have curated a set of benchmark cases through collaborative effort with several research groups [16].The goal is to first verify the physicochemical models and the numerical schemes.Once confidence is established over these fundamental aspects, the accuracy and limitations of the IB methods is addressed.We have selected setups that are sufficiently challenging for the methods under assessment, and simple enough to be readily reproduced by others to incentivize collaboration.For IB methods on Cartesian grids, curved geometries were selected to include the entire angular range of fluid-solid interfaces in two dimensions.These cases include strong thermal gradients near cold isothermal walls as well as gas-surface interactions such as catalytic reactions and ablative surface blowing.
The benchmark cases are summarized in Table 1.The first four cases serve as the verification of the implementations for chemistry, transport properties, the catalytic boundary conditions, and the numerical schemes for shock capturing.Established validation experiments are chosen as the final benchmark cases: the fifth, sixth and seventh cases are 2-D cylinder flows of Knight et al. [3] with inert adiabatic, inert isothermal, and catalytic isothermal surfaces.As the eighth benchmark, we discuss results for an ablative TPS sample geometry under plasma wind tunnel conditions, for which reference experimental data is provided by Helber et al. [53].2, the system is left free to time-march towards the equilibrium state according to chemical mechanisms from Park [54,55].The solutions provided by all three solvers are shown in Fig. 2. Dissociation of N 2 , O 2 and the resulting formation of NO, N, and O can be seen.The code-to-code agreement is excellent.Mutation ++ , and are exactly equal for all solvers.Therefore, mainly the differences in the implementation of the driving force and boundary conditions are assessed.The setup consists of a 1-D tube with isothermal end walls at different temperatures.The initial and boundary conditions are provided in Table 3.The mixture composition and reaction mechanisms are the same as in the 0-D reactor case.The tube is 3 mm long.It should be pointed out that the computational meshes in US3D and INCA solvers have 100 cells, whereas CHESS results [16] used 400 cells.It has been verified that the US3D and INCA solutions are grid converged on the mesh with 100 cells.
0.02 1000 800 4800 0.0 0.767 0.233 In this test case the temperature gradient leads to chemical reactions, which in turn drive mass diffusion.Temperature and mass fraction distributions along the tube are presented in Fig. 3. INCA results have been obtained by both Fick's law and Stefan-Maxwell diffusion models.However, for this test case, differences seem to be negligible between the two.Overall, US3D results are matched perfectly with INCA, while slight differences are observed for the mass fraction distributions predicted by CHESS, even though the temperature profiles match exactly.

1-D Catalytic Diffusion Problem
This test case verifies the catalytic boundary condition implementation for a simple [N 2 , N] binary mixture along a 1-D tube, for which an analytical solution exists and is derived in the appendix.Setup conditions are given in Table 4.The length of the tube is 0.2 m.One side of the tube at x = 0.0 m is at reservoir conditions, while at the other, at x = 0.2 m, a catalytic wall boundary condition is imposed.The catalytic wall promotes the recombination of nitrogen, that is, N + N → N 2 .The reaction rate is controlled by the recombination coefficient γ through Eq. 9.  Results obtained with US3D, INCA, and CHESS [16] are compared with the analytical reference solution in Fig. 4. Naturally, for higher values of the recombination coefficient γ, mass fraction of molecular nitrogen at the wall increases, and reaches unity for the fully catalytic case with γ = 1.0.All numerical predictions are in excellent agreement with the analytical solution.The previously noted difference for the CHESS solver in the diffusion problem is not observed here as the diffusion of species are driven predominantly by the surface reactions.

1-D Shocktube
The Riemann problem of Grossman and Cinnella [56] is used to evaluate shock capturing.The unit domain, x = [0, 1] m, is spatially discretized by 600 cells, in line with the reference resolution.The diaphragm separating the two initial states is set at the midpoint of the tube.The initial conditions for the two states are given in Table 5.Air with 5-species is initially considered at thermodynamic equilibrium.The reaction mechanism is taken from an earlier work of Park [57], to match with the reference [56].Grossmann and Cinnella applied a thermal nonequilibrium model; however, we have performed tests with Park's two-temperature model [19] and found no significant differences between the translational and vibrational energy modes.Therefore, we show results that have been obtained with a thermal equilibrium assumption.
0.0 9000 195256 0.0 300 10000 Fig. 5 shows pressure and density profiles 99 µs after the initial state.Mass fractions are given in Fig. 6.The contact discontinuity and the shock wave traveling in the positive x direction as well as the expansion traveling in the opposite direction are well captured.The peak in density after the shock also matches perfectly with the reference results without any oscillations.Predictions of US3D and INCA for the mass fractions are also in excellent 'LVWDQFH>P@ 3UHVVXUH>N3D@ 86' ,1&$ &+(66,% *URVVPDQ &LQQHOOD agreement with the reference results of Grossmann and Cinnella, see Fig. 6a.The minor differences between the solvers in their sharp representation of the discontinuity is shown in the close-up view in Fig. 6b.CHESS [16] predicts a slightly higher N 2 mass fraction, and accordingly less atomic nitrogen, than US3D and INCA.

2-D Cylinder
2-D cylinder flows are used for the validation of surface heat flux calculations under inert and catalytic wall conditions.Knight et al. [3] have presented an assessment of five different CFD codes from participating institutions with respect to reference experiments conducted at the high-enthalpy shock tunnel of the German Aerospace Center (DLR) [58].The experiment investigates the flow past a cylinder with a radius of 0.045 m exposed to a reported total enthalpy of 22.4 MJ/kg.The experimental setup is numerically replicated by imposing the inflow conditions given in Table 6 on the left boundary.Symmetry is imposed along the stagnation line, and the outer boundaries are set as non-reflecting outlets.The reaction mechanism employed for the 5-species air model is taken from Park [54,55].As remarked by Knight et al. [3], there appears to be a large variation in the results from different solvers, especially regarding the treatment of the surface.To study this sensitivity, three different surface conditions are tested in the following sections: two inert cases with adiabatic and isothermal conditions, and a third case with a fully catalytic isothermal wall.'LVWDQFH>P@ 0DVV)UDFWLRQV>@ 1 2 In the following assessment of IB methods, "BF" is used to denote reference results obtained on body-fitted grids with US3D, "IB-CC" is used for the cut-cell IB method of INCA, and "IB-GP" refers to the ghost-cell IB method of CHESS [16].

Inert Adiabatic Wall
The temperature and species mass fractions along the stagnation line are presented for the adiabatic case in Fig. 7. Shock stand-off distance and the dissociation of molecular nitrogen and oxygen in the shock layer are predicted in very good agreement between all methods.The fundamental differences in the implementation of the adiabatic wall boundary condition have no noticeable effect on the results.This is in line with the expectation that truncation and conservation errors are small in the absence of strong gradients.'LVWDQFHIURPWKHVXUIDFH>PP@ 7HPSHUDWXUH>.@%) ,%&& ,%*3

Inert Isothermal Wall
For the same inflow conditions, an isothermal wall boundary condition with a wall temperature of 300 K is imposed on the cylinder surface, in accordance with the specifications by Knight et al. [3].The numerical predictions for the stagnation line temperature and mass fraction distributions are plotted in Fig. 8. Results obtained with the BF and the IB-CC methods match almost exactly, including the steep temperature and species variations in the boundary layer.Results obtained with the IB-GP method, on the other hand, show a significant difference in the shock stand-off distance.This could be attributed to the non-conservative formulation of the ghost-cell IB methodology.Mass conservation errors could manifest as an unphysical blowing from the surface.Consequently, the shock stand-off distance is increased and the whole flow field is modified.The adiabatic case is less affected by these conservation errors because it has much smaller temperature and density gradients near the wall.The IB-CC method handles large temperature and density gradients at isothermal walls much better, because it uses a strictly conservative IB method.In Fig. 9, surface pressure and heat flux distributions are compared with the experimental measurements from Knight et al. [3] and also with the numerical simulations of Nompelis from the same publication.All methods accurately predict the pressure distribution.Heat flux predictions of the BF and the IB-CC methods are in very good agreement.They match the experimental measurements better than the numerical simulations of Nompelis.Slight differences in heat fluxes are expected to be due to the differences in grid resolutions at the surface.Grid convergence studies have been carried out with both the BF and the IB-CC methods as summarized in Table 7 and showcased for the variation in heat fluxes in Fig. 10.Four levels of resolution are considered with the minimum cell size at the surface approximately halving with each step.For both solvers, results obtained on the medium-fine resolution mesh are considered grid converged, as they are essentially identical to the results obtained on the fine meshes.For these grids, smallest cell size near the wall is 1.0 × 10 −7 m for the BF method and 6.25 × 10 −7 m for the IB-CC method.An interesting observation is that IB-CC method under-predicts the heat flux on coarse meshes, as intuitively expected, whereas the BF method over-predicts the heat flux on coarse meshes.This difference in the convergence trend is a sign of complex interactions between transport and chemistry.
Table 7: Grid resolution and stagnation point details for the inert isothermal 2-D cylinder case, where ∆h w is the effective wall-normal cell size at the wall, p 0 is the stagnation point pressure, and q w,0 is the stagnation point wall heat flux.For this case, the IB-GP method is not able to predict the heat flux correctly.A similar underprediction has also been reported in literature for another ghost-cell IB-AMR solver by Brahmachary et al. [12], where the issue has been linked to the reconstruction of temperature by linear interpolation.However, the cut-cell IB method also resorts to second-order reconstruction schemes and can predict the heat flux correctly.Therefore, we attribute the observed deficiencies to conservation errors incurred through the ghost-cell IB method.This hypothesis is further discussed in Section 5, where results obtained with two independently developed ghost-cell IB methods are presented.

Catalytic Isothermal Wall
Exothermic catalytic reactions enhance the surface heat flux through a diffusive heat flux contribution.A fully catalytic wall (γ = 1.0) at the same temperature of 300 K is considered, as Karl et al. [58] state that this boundary condition is closest to what they have assumed for the experiments.It is, however, reasonable to assume that fully catalytic conditions were only achieved for a short duration at the beginning of the experiment.The fully catalytic boundary condition imposes the recombination of all atoms impinging on the surface, while still respecting the physical limits set by species diffusion.
The species mass fractions along the stagnation line and the total surface heat flux distributions are shown in Fig. 11.Predictions of the BF and the IB-CC methods are in excellent agreement, both in terms of species concentrations and surface heat fluxes.Because the cold wall itself already promotes recombination reactions in the boundary layer, accounting for catalysis leads only to a minor increase in the heat fluxes, which remain within the experimental uncertainties.It is therefore difficult to draw conclusions on the effective value of the recombination coefficient in the experiment.
Another interesting observation could be made by comparing the level of agreement between the BF and the IB-CC results for the heat fluxes at the inert wall shown in Fig. 9 and with the fully catalytic wall shown in Fig. 11.Taking the BF method as reference, it is seen that at the stagnation point, results of the IB-CC method with the inert wall are 2.7% lower, while with the fully catalytic wall the difference is only 0.2%.It can be argued that this better agreement is mostly associated with the dominant nature of the catalytic boundary condition.Our previous comments regarding the differences seen in the IB-GP method for the inert isothermal case apply here as well.
To complete the analysis, contour plots are presented for Mach numbers in Fig. 12, for temperatures in Fig. 13, and for atomic nitrogen concentrations in Fig. 14.These contour plots further confirm the preceding quantitative discussions by once again reflecting the excellent agreement between the BF and the IB-CC methods.From the trace of the sonic line, to the peak shock temperatures, and to the extent of the species boundary layer marked by nitrogen accumulation, the results are in perfect agreement.

2-D Ablator
To validate the ablative boundary condition implementation and to assess the IB methods for GSI, a subsonic plasma wind tunnel experiment conducted at the von Karman Institute for Fluid Dynamics (VKI) by Helber et al. [53] is considered.The experiment exposes a graphite sample with a hemispherical nose of radius 25 mm and a downstream extension of 250 mm to nitrogen plasma.The sample undergoes ablation through nitridation reactions The nitridation efficiency was calibrated based on these particular experiments [53].The simulations discussed in the following include mass blowing due to ablation, but do not account for the very slow shape change of the sample.First, we reproduce the experiment numerically using the BF method.For these simulations, a 9-species nitrogen-carbon mixture is considered, including free electrons and ionized species.These simulations yielded a stagnation point mass blowing rate of 3.41 g/m 2 s, which is within the experimental uncertainty range set by 2.8864 ∓ 0.965 g/m 2 s.This validates the ablation model based on Eq. 13.
Having confidence in the ablation model and its implementation in the BF method, the experimental test case is simplified to a 2-D geometry without ionized species to reduce the computational cost and to avoid straying too far from the objective of evaluating immersed boundary methods for    an ablative boundary condition.Freestream conditions of this 2-D case are given in Table 8.A 6-species mixture of [N 2 , N, CN, C 3 , C 2 , C] is considered with chemical mechanisms from Olynick et al. [59].For all methods, the grid resolution at the wall is 1 × 10 −5 m in the wall-normal direction.u The mass fractions along the stagnation line and the mass blowing rates over the wall are shown in Fig. 15.Mass fractions for C 3 are not seen as they are almost zero.Predictions of the BF and the IB-CC methods agree well with each other.Overall, the production of CN at the wall and the dissociation of it through gas-phase reactions to form atomic nitrogen are well captured.Mass blowing rates from the BF and the IB-CC methods are also in very good agreement.The IB-GP method show noticeable discrepancies for the mass fractions along the stagnation line and for the surface mass blowing rates.Despite the apparent quantitative mismatch, also the IB-GP method captures the profiles qualitatively well in the absence of strong gradients near the wall.
Temperature and atomic nitrogen contours for the simulations with the BF and the IB-CC methods are shown in Figs.16 and 17. Results of both methods agree very well on the thermal gradient over the surface and on the recombination of nitrogen as temperature drops.

On the Importance of Conservative Boundary Conditions
In the previous section, we have demonstrated that INCA with its cut-cell IB method on block-Cartesian AMR meshes performs on par with the reference solver US3D employing body-fitted meshes.The ghost-cell IB method of CHESS yields similar accuracy for cases with adiabatic walls; however, it cannot predict the heat flux at strongly cooled walls.We have attributed these inaccuracies to mass conservation errors as this is the most striking difference between ghost-cell methods and the strictly conservative cut-cell approach.However, the three solvers clearly differ also in several other aspects, such as the numerical schemes used for advection and diffusion driving forces.It is therefore unclear whether the observed deficiencies are inherent to the ghost-cell method or a particular implementation.To further corroborate the superiority of a conservative cut-cell (or cut-element) IB methodology, we have also applied the ghost-cell method of INCA for selected cases.By switching off the special flux treatment in cut-cells [51] employed in the preceding sections, a standard ghost-cell method is obtained that only relies on the extrapolated fluid solutions near the boundary as described in Section 3. Therefore, mass, momentum, and energy conservation are not exactly satisfied.This ghost-cell method has nominally the same order of convergence as the conservative cut-element method of INCA.The comparison of the two methods is shown for the two most challenging benchmark cases in Fig. 18, where the previous cut-cell based results from the INCA solver are denoted by "INCA-CC" and the ghost-cell based results are denoted by "INCA-GP".Similarly, results with the ghost-cell IB method of the CHESS solver are denoted by "CHESS-GP".It can be seen that regardless of the various differences between INCA and CHESS, in both cases the ghost-cell methods are unsuccessful in predicting the surface heat fluxes and the mass blowing rates.For the 2-D cylinder case with an isothermal wall, the heat flux prediction of the INCA-GP method is closer in magnitude to the INCA-CC method and to the body-fitted reference from US3D than to the results obtained with the CHESS ghost-cell IB method; however, both ghost-cell methods give clearly wrong results.For the ablative case, the ghost-cell based method of INCA yields a very similar overprediction as the IB method of CHESS.ents of the conservative variables are very large.This explains why errors manifest at cold walls and not at adiabatic walls.

Conclusion
We have evaluated the accuracy of immersed boundary methods for atmospheric entry conditions, including the influence of chemical nonequilibrium and gas-surface interactions.The benchmark cases have considered the accurate modelling of gas chemistry, mass diffusion, surface catalysis, and surface mass blowing due to ablation.
Computational results obtained with the cut-element IB method in the AMR solver INCA are in an almost perfect agreement with the reference data for all considered cases, and as accurate as the results obtained with US3D on body-fitted meshes.Particularly for surface heat flux and mass blowing rate predictions, the benefit of an IB method that strictly conserves mass, momentum, and energy, such as the cut-element method in INCA, is clearly demonstrated in this study.After comparing this method with two non-conservative ghost-cell methods implemented in INCA and an independently developed solver, we saliently remark that numerical anomalies causing mispredictions of sensitive surface quantities can occur when using non-conservative IB formulations.
CFD solvers that provide automatic mesh generation and adaptation to represent detailed and moving geometries with IB methods have many promising advantages, but the accuracy of the numerical schemes used for predicting surface quantities must be analyzed rigorously before they can be used for predictive simulations.The selection of a set of well-defined test cases by mutual collaboration between research groups is crucial in converging to a robust consensus on the prediction of these surface states.To that end, this paper establishes such a set of fundamental benchmark cases with reacting surfaces, which can be used for the verification and validation of hypersonic flow solvers, while assessing the accuracy of immersed boundary methods for atmospheric entry applications.with η as the spatial coordinate.Solving for x N 2 yields with C 1 and C 2 as integration constants to be found through the boundary conditions.Firstly, by knowing that (X N 2 ) η=0 = 0 at the free-stream reservoir Secondly, by equating the diffusion flux to the chemical production rate at the wall, (J N 2 = ωN 2 ) η=L , which gives where k B is the Boltzmann constant.The last expression can be solved iteratively through the Newton-Raphson method.The solution describes the species distribution as a function of spatial variable η.

4. 1
. 0-D Reactor The first study verifies the chemical source term implementation by considering 5-species air, [N 2 , O 2 , NO, N, O], in an adiabatic reactor.Starting from the chemical nonequilibrium (CNE) initialization provided in Table

Figure 2 :
Figure 2: Evolution of mass fractions for 5-species air in the 0-D reactor case.

Figure 3 :
Figure 3: Comparison of (a) temperature and (b) mass fraction distributions for the 1-D diffusion case.

Figure 4 :
Figure 4: N 2 mass fractions for different recombination coefficients γ for the 1-D catalytic diffusion problem.

Figure 5 :
Figure 5: Comparison of (a) pressure and (b) density distributions for the 1-D shocktube case of Grossman and Cinnella [56].

Figure 6 :
Figure 6: Comparison of mass fraction distributions for the 1-D shocktube case.

Figure 7 :Figure 8 :
Figure 7: Comparison of (a) temperature and (b) mass fractions along the stagnation line for the inert adiabatic 2-D cylinder case.

Figure 9 :
Figure 9: Comparison of surface (a) pressures and (b) heat fluxes for the inert isothermal 2-D cylinder case by Knight et al. [3].

Figure 10 :Figure 11 :
Figure 10: Grid independence studies with (a) the BF and (b) the IB-CC methods considering surface heat fluxes for the inert isothermal 2-D cylinder case by Knight et al. [3].

Figure 12 :
Figure 12: Comparison of Mach number contours for the catalytic isothermal 2-D cylinder case by Knight et al. [3] obtained with (a) the BF and (b) the IB-CC methods.The sonic line is indicated with the bright yellow line.

Figure 13 :
Figure 13: Comparison of temperature contours for the catalytic isothermal 2-D cylinder case by Knight et al. [3] obtained with (a) the BF and (b) the IB-CC methods.

Figure 14 :
Figure 14: Comparison of atomic nitrogen mass fraction contours for the catalytic isothermal 2-D cylinder case by Knight et al. [3] obtained with (a) the BF and (b) the IB-CC methods.

Figure 15 :
Figure 15: Comparison of (a) mass fractions along the stagnation line and (b) mass blowing rates over the surface for the 2-D ablator case.

Figure 16 :
Figure 16: Comparison of the temperature contours for the 2-D ablator case obtained with (a) the BF and (b) the IB-CC methods.

Figure 17 : 2 @Figure 18 :
Figure 17: Comparison of the atomic nitrogen mass fraction contours for the 2-D ablator case obtained with (a) the BF and (b) the IB-CC methods.

Table 1 :
Summary of studied cases.

Table 3 :
Setup conditions for the 1-D diffusion case.

Table 4 :
Setup conditions for the 1-D catalytic diffusion case.

Table 5 :
Initial conditions for the 1-D shocktube case.

Table 6 :
Freestream conditions for the 2-D cylinder case.