The variable‐extended immersed boundary method for compressible gaseous reactive flows past solid bodies

An immersed boundary method (IBM) has been developed to handle the solid body embedded flowfield simulation for compressible reactive flows, paving the way of application for a wide range of fluid–solid interaction problems. Previously, the Brinkman penalization method (BPM), originated from porous media flows, has been successfully used for incompressible Navier–Stokes equations by adding penalization terms to momentum equations. However, it is nontrivial to solve the compressible form due to the penalized continuity equation that usually poses severe numerical stiffness. In order to circumvent this issue, an extending procedure for relevant variables from the fluid to solid domain is considered, by analyzing the ordinary differential equations remained after operator splitting. Density can be then determined with the help of an equation of state (EoS). Meanwhile, efforts of enforcing the Neumann boundary condition, for example, the adiabatic wall condition, on the fluid–solid interface can be minimized by extending temperature across the interface directly. One more advantage of the extending step lies in that it can quickly reach a steady state when performed within a narrow band around the interface. Implemented into an adaptive Cartesian grid‐based flow solver for compressible Navier–Stokes equations with chemical reaction source terms, the present variable‐extended IBM is validated by numerical examples ranging from single‐species nonreactive to multispecies detonative flows in one‐ and two‐dimensional domains. Numerical results show (1) the successful specification of slip or nonslip, adiabatic or isothermal wall condition on the fluid–solid interface and (2) loss of total energy in the original BPM being avoided and the numerical accuracy being improved especially for energy‐sensitive reactive flows.

complex geometries and the associated boundary conditions. Two main categories of methods are usually considered to handle complex geometries: body-fitted grid-based methods and immersed boundary methods (IBMs). Body-fitted methods rely on structured or unstructured grids that are conformal to the complex boundaries. To generate a mesh of good quality, one should be very careful even for some simple geometries, which leads to the mesh generation process being prohibitively expensive. Besides, for a given spatial discretization scheme, the order of accuracy on structured or unstructured grids is always lower than on Cartesian grids. 1 Alternatively, one can choose to perform simulations on simple Cartesian grids and to impose immersed boundary conditions upon the fluid. Two essential advantages of this idea are its easy implementation and relatively straightforward extension from stationary body to more general moving body cases without mesh regeneration.
A large group of IBMs belong to the direct-forcing IBM. It on one hand mimics the behavior of body-fitted methods that use ghost cells, and imposes the boundary conditions directly on the immersed boundaries. [2][3][4][5][6][7][8] Multidimensional interpolation or extrapolation in the vicinity of the fluid-solid interface is usually needed in order to determine the image-point state in various cases and thus the ghost-cell state, making use of intersection points on the interface to impose exact boundary conditions. 9 On the other hand, some cut-cell methods [10][11][12][13] confront the obstacles about the specification of the type or local geometry of a cut cell (especially in high-dimensional situations) as well as the fact that small cut cells might lead to accuracy losses and/or strict timestep restrictions. 14 As one another type of IBMs, the Brinkman penalization method (BPM) solves the fluid and the solid domain simultaneously by introducing fictitious (penalization) terms in a unified system of governing equations. In this way, there is no need to impose the interface conditions explicitly since they are automatically solved from the governing equations. The BPM originated from the idea of considering the solid body as a porous medium with a very low permeability and the fluid as a medium with an infinite permeability by Arquis and Caltagirone, 15 and became available for incompressible Navier-Stokes equations first by Angot et al. 16 Liu and Vasilyev 17 started to extend the BPM to the regime of compressible viscous flows by penalizing the continuity equation in addition to the momentum and energy equations. They also pointed out that only penalizing the momentum and energy equations can produce unsatisfactory results, mostly due to nonphysical wave transmissions into solid bodies, resulting in considerable energy and mass losses in reflected waves. Boiron et al. 14 further considered the BPM for large Mach number flows. Anand et al. 18 investigated the BPM in combination with a high-order discontinuous Galerkin scheme, and stated that the modeling error (e.g., ( 1∕2 ) for resolved boundary layers 17 ) can be minimized with sufficiently small permeability when holding the porosity = 1. Piquet et al. 1 utilized an operator-splitting scheme to separately treat the penalization terms in a semiimplicit manner. All the works above focused the discussion on the isothermal wall condition or Dirichlet boundary condition for temperature or total energy. Besides, to apply the Neumann boundary condition for temperature, Browndymkoski et al. 19 proposed a characteristic-based volume penalization (CBVP) method that takes into account adiabatic walls and mixed boundary conditions. With CBVP, Hardy et al. 20 simulated weakly compressible reacting gas-particle flows using a predictor-corrector scheme on a staggered mesh. In a general sense, severe numerical stiffness in the penalized continuity equation and nontrivial imposition of the Neumann boundary condition are two main bottlenecks in BPM. 9,19 Aiming at the simulation of fully compressible gaseous reactive flows past solid bodies, the present study begins with the penalized compressible Navier-Stokes equations with chemical reaction source terms, including the penalization of the continuity equation. To overcome the numerical stiffness resulted from the penalized continuity equation, an extending step for relevant primitive variables from the fluid to solid domain is considered to replace solving the ordinary differential equations (ODEs) of the penalization terms after operator splitting. Density inside the solid body can be then determined with the help of an EoS. Meanwhile, by extending temperature across the interface directly, which is analogous to the treatment in the body-fitted methods, the Neumann boundary condition (e.g., the adiabatic wall condition) can be applied in a considerably simple and robust manner. Special attention to imposing the velocity condition on the immersed boundary is also paid to avoid the loss of total energy by the original BPM. To minimize the computational cost and numerical complexity, the present extending operation on one hand is only performed within a narrow band around the fluid-solid interface in the normal direction from the fluid to solid domain; it on the other hand can converge to a smooth steady-state solution unconditionally, given the nature of the advection equation (i.e., the extending equation) with unit velocity along the normal direction. Having implemented the present method into an adaptive Cartesian grid-based flow solver for compressible reactive Navier-Stokes equations, we carry out a variety of numerical examples from single-species nonreactive to multispecies detonative flows in one-and two-dimensional (1D and 2D) cases.
The article is organized as follows. In Section 2, we first introduce the fully compressible Navier-Stokes equations with chemical reaction source terms as well as boundary conditions on the fluid-solid interface. Then the BPM is considered to enforce such boundary conditions into the compressible reactive Navier-Stokes equations. The present variable-extended IBM is followed. In addition, we specify the spatial and temporal discretization schemes in our adaptive Cartesian grid-based flow solver, and only consider one-way fluid-solid interaction in the scope of this study since the extension to moving body dynamics is straightforward. In Section 3, we examine the present method by extensive numerical examples. Conclusions are drawn in the last section.

Governing equations
Let Ω be the computational domain containing N fixed regular obstacles n s , 1 ≤ n ≤ N, and we set Here, Ω s is the closed region occupied by the solid bodies and Ω f is the fluid domain. For the fluid domain Ω f , we consider the fully compressible Navier-Stokes equations coupled with gaseous finite-rate nonequilibrium chemistry, together with appropriate boundary conditions on the solid body surface (or fluid-solid interface) n s and on the boundary of the computational domain Ω. For simplicity and without loss of generality, the system of equations in 2D space is considered and has the conservative form where , and are the vector of conserved variables, the inviscid flux vectors, the viscous flux vectors in the x and y directions and chemical reaction source terms, respectively. Here , u = [u, v] T , p, T, and e t denote the density, velocity, pressure, temperature, and the specific total energy, respectively, and e t = e + 1 2 (u 2 + v 2 ) including the specific internal energy e. y i is the mass fraction of the ith species in the ns-species reactive gas mixture, and constrained by mass positivity and conservation, that is, y i ≥ 0 and ∑ ns i=1 y i = 1. In the viscous flux vectors, ij is the shear stress defined by with being the dynamic viscosity of the fluid. Viscosity can be calculated using Sutherland's law for single species in combination with a proper mixing rule such as the Wilke's semiempirical formula. 21 The heat flux is where k is the thermal conductivity and molecular diffusion terms can be approximated by using with Sc = 0.5 being the Schmidt number. In the source terms,̇i represents rate of change of the ith species concentration due to nr finite-rate chemical reactions, that is, and we havėi with l = y l . When the flow is inert or nonreactive without activating chemical reactions, source terms are replaced by a zero vector. To close the system, an EoS of the form is used for the ideal gas mixture, with W i denoting the molecular weight of the ith species and R u the universal gas constant.
On the surface of each solid body n s , the fluid velocity should satisfy the no-slip condition given the moving solid body velocity u n s . u n s = 0 corresponds to a stationary solid body. For variables such as density, pressure, temperature, species mass fractions, and so forth, the Dirichlet boundary condition as (given the value of any variable on the solid body surface), or the Neumann boundary condition as is usually imposed. For example, isothermal and adiabatic wall boundary condition for temperature on the solid body surface correspond to these two cases, respectively, with = 0 in the latter one.

Brinkman penalization method
One very easy-to-implement method to enforce the nonslip velocity condition on solid body surfaces in IBMs is BPM, which considers the solid body as a porous medium with a very low permeability and the fluid as a medium with an infinite permeability. It was originally developed for simulating incompressible flows 16 and imposing the Dirichlet boundary condition for other variables in addition to the nonslip condition for velocity in Equation (10). For the fully compressible Navier-Stokes equations, Brinkman penalization terms are added, giving and where is the porosity with 0 < ≪ 1, is the permeability with 0 < ≪ 1 and is the characteristic (or mask) function of the solid domain Ω s . It is readily to see that the penalization terms are effective only in the solid domain with = 1. It is worth noticing that the above system of equations is solved for both the fluid and the solid domain. Using the mask function , the two domains are automatically defined and treated in a smooth and unified manner. This is a merit of BPM since it inherently skipped the geometric details of the fluid-solid interface inside a cut cell and avoids dealing with the various types of cut cells in high-dimensional situations with other IBMs. Extensive studies 17,22 have shown the controllable convergence of the penalization method by prior parameters such as by decreasing and , or increasing the interface grid resolution. Based on operator splitting, the penalization terms (P) can be separated from the convection-diffusion reaction terms (C). 1 For example, considering the second-order Strang splitting, 23 we have the following fractional steps as Therefore, in terms of the penalization step, a series of ODEs are left, for example, for the momentum and energy equations, and an analytical solution exists for such a type of ODE 1 as for a time interval Δt = t n+1 − t n , in spite of that might be a very small number. The penalized continuity equation, however, takes a more complicated form if ≠ 1 according to Equation (14) as Since 0 < ≪ 1 leads to 1 ≫ 1, the penalized continuity equation might exhibit very severe numerical stiffness in terms of a finite timestep of Δt. In such a case, a very tiny local timestep on the solid body surface is necessary to overcome the stiffness, which inevitably limits the size of the global timestep and leads to the low overall computational efficiency. Therefore, some studies 14,18 simply eliminate the penalization terms in the continuity equation by assuming = 1 at the cost of losing some numerical accuracy. It is also possible to use implicit schemes to deal with the numerical stiffness while the increased computational complexity must be compromised with the decreasing size of . In some practice, the CFL number cannot yet be greater than 1.0 even though the implicit scheme such as LU-SGS 24 is employed for reactive flows. The severe numerical stiffness in the penalized continuity equation is one disadvantage of BPM. 9 Another one lies in the nontrivial application of the Neumann boundary condition in Equation (12). 19

Variable-extended IBM
With the help of Equation (9), we see that the penalized density can be obtained directly through the EoS using other primitive variables such as pressure, temperature and mass fractions, instead of solving the stiff penalized continuity equation in Equation (19). In Reference 19, the penalization equation for temperature with the Neumann boundary condition (e.g., = 0 in Equation (12)) in the solid domain ( = 1) is given as where n = [n x , n y ] T is the unit normal vector of the solid body surface pointing outward to the fluid domain. We notice that Equation (20) is similar with the extending equation in ghost fluid methods 25,26 that has the form of an advection equation with being a pseudo time here. It is readily to see that a steady state, corresponding to T n = 0 on the solid body surface and inside the solid domain, can be derived from the extending equation and is one solution to Equation (20), regardless of how small the size of is or how relatively large the size of Δt is. More extending equations are available for pressure and mass fractions if the Neumann boundary condition is considered. Note that the direction of extending should be pointed from the fluid domain towards the solid domain.
Remark 1. The extending procedure resembles the ghost-cell methods that need to populate the ghost-cell states by reflecting the interpolated image-point states inside the fluid domain in combination with the exact boundary condition of the intersection points on the fluid-solid interface. The numerical complexity results from the interpolation (or extrapolation) schemes and the treatment of different types of cut cells along with intersection points. However, the present extending operation derives from the penalization ODE after operator splitting, and it can converge to a smooth steady-state solution unconditionally. It therefore spares the efforts of considering the geometric complexity of local interface especially in the existence of cut cells, image points and intersection points, and inherits such a critical merit of BPM.
Solving the extending equation such as Equation (21) is computationally inexpensive with respect to the narrow band of ghost cells on the solid body surface and inside the solid domain, as depicted in Figure 1. The width of the narrow band mainly depends on the length of the interpolation stencil for high-order reconstruction schemes in the calculation of inviscid fluxes. For example, given the ghost cell (i, j) which is in direct contact with two fluid cells and the local normal vector n in the figure, one can easily find the extending relies on the neighboring cells (i − 1, j) and (i, j + 1). In general, considering a linear approximation, we have the explicit expression where For the remaining cells not in direct contact with fluid cells, using Equation (22) in an iterative manner will result in the desired steady state quickly within the narrow band. The computational cost of such an extending procedure over the narrow-band area only is much less than solving the penalized Navier-Stokes equations over the entire solid domain in the original penalization method.
Remark 2. The iterative algorithm here is robust and can maintain the smoothness of the flowfield extended from the fluid domain to the solid domain around the interface. The resulting accuracy is easily adjustable by using different schemes to approximate the gradient term in Equation (21).
Similar with extending temperature from the fluid domain to the solid domain, pressure and mass fractions of species can be also extended subsequently. The EoS will then supplement other variables needed in the ghost cells and a complete set of data used for the inviscid/viscous flux calculation can be constituted.
One more concern in the compressible reactive flow simulation exists, that is, even though the penalized velocity inside the solid body will quickly recover the solid body velocity as in Equation (10) the total energy. The loss of energy does not matter much for incompressible or weakly compressible flows. However, for compressible flows, especially for the energy-sensitive reactive flows, such a loss possibly gives rise to inaccurate time-dependent solutions, for example, a slowly propagating reacting front which will be given a description in a following numerical example. Using Figure 2(A) and assuming the solid wall is stationary (u w = 0), we can find that while the interior cell in the fluid domain (i, j) has a positive velocity of u i > 0, its image ghost cell (i + 1, j) inside the solid wall ought to have a velocity of u g = 0 by the penalization method. However, we can recall that in the body-fitted grid-based method, an exact nonslip wall condition essentially gives u g = −u i such that the total energy e t = e + 1 2 u 2 on both sides of the wall equals each other. In the penalization method, an energy loss of 1 2 u 2 is a by-product, which is likely to contaminate the accuracy of the compressible reactive flow simulation. Therefore, we continue to replace the penalization step by an extending step for velocity, similar with that for temperature in Equation (22), to have u n+1 i+1,j in Figure 2(A), for example. This is followed by an additional step of changing the sign of the postextending velocity, that is, to fulfill the nonslip condition. In some cases, the slip velocity condition for inviscid flows is also in consideration, then we obtain the image velocity vector with respect to the interface using Regarding the above extending procedures, we emphasize that the present method merely requires (1) the shortest distance from each solid body cell to the fluid-solid interface to determine the narrow band of ghost cells and (2) the unit normal vector of each ghost cell to the fluid-solid interface. As the fluid-solid interface is usually defined by the mask function , a levelset function is widely used to play such a role. However, the levelset function involves a costly on-the-fly reinitialization procedure when the fluid-solid interface is moving. In this work, we can simply use a series of Lagrangian points (markers) to outline the closed fluid-solid interface in an arbitrary user-defined polygon shape, 27 so that geometry-related information including the shortest distance and the interface normal can be conveniently provided in an efficient and robust manner. As shown in Figure 2(B), the Lagrangian points are independent of the computational cells/points in the background grid and queued in the counter-clockwise direction. Given the ghost cell G inside the solid domain, the shortest distance from G to the interface is the shortest one of distances from G to edges P k P k + 1 using a simple sorting algorithm. Note that an eligible distance candidate should satisfy where np is the number of points/vertices of the polygon, and P k + 1 ← P 1 when k = np. Besides, perpendicular to the nearest edge is the normal vector from cell G we desire. Other properties such as the centroid of the solid body, total mass, and inertial moment with respect to the centroid are all available if needed.

Spatial and temporal discretization in fluid domain
Discretizing the compressible Navier-Stokes equations in Equation (2) on an adaptive Cartesian grid, 25 flow solver in the present study can conveniently employ high-order low-dissipation shock-capturing schemes to reconstruct the inviscid convective flux terms, as well as high-order central difference schemes to compute the viscous diffusive fluxes. In this study, the fifth-order WENO-LLF 28 finite difference scheme for multispecies reactive flows, based on characteristic decomposition 29 and upwind flux splitting, and a simple fourth-order central difference scheme are used, respectively. To balance the overall accuracy in time and the computational cost, the temporal integration utilizes the strong stability preserving second-order explicit Runge-Kutta scheme for a timestep Δt constrained by the CFL condition.
Remark 3. Despite the present variable-extended IBM being more like a smooth interface method, its capability of sharp interface capturing can be alternatively strengthened by improving the spatial resolution near the interface with adaptive grid refinement.

Fluid-solid interaction
In the scope of this article, we simply consider the stationary solid body embedded in the computational domain. Hence, the fluid-solid interaction is a one-way coupling here. That is to say the solid body prevents the flow from passing through directly but is not influenced by the fluid. However, further implementation of moving solid bodies and the related dynamics is straightforward. For example, given the flowfield information around the solid body, fluid-solid interaction can be considered by imposing forces such as pressure and viscous forces from the fluid upon the solid body surface. Each vertex of the solid body polygon can be used as a sampling point that measures the local pressure and viscous forces. The integrated forces and torques will then be used to determine the motion of rigid bodies if the rigid-body dynamics is considered, see more details in References 30,31.

NUMERICAL RESULTS AND DISCUSSION
In this section, we attempt to validate the present variable-extended IBM using several numerical examples concerning simple single species (ideal air) or multispecies reactive flows with finite-rate nonequilibrium chemistry past the solid body obstacle in 1D or 2D space. For solving the ODE system of chemical kinetics, an efficient and robust explicit ODE solver, CHEMEQ2, 32 is employed. The kinetic mechanism which describes the species and reactions will be specified in relevant examples.

Case 1: One-dimensional acoustic wave reflection
To assess the effectiveness of the present variable-extended IBM in capturing the reflective nature of a solid wall, we use the reflection of a 1D acoustic wave at the interface between the fluid and a solid obstacle. The entire computational domain Ω is [0,0.6]. The fluid occupies Ω f = [0, 0.5] and the solid occupies Ω s = [0.5, 0.6] with the interface being located at x = 0.5. The initial conditions include the density, velocity and pressure perturbations of the Gaussian distribution where the wave amplitude is = 10 −3 , to be applied upon the spatially uniform field of = 1, u = 0, p = 1 , and c = 1 being the speed of sound. The fluid material is air with = 1.4. Theoretically, wave reflection at the wall x = 0.5 should be perfectly symmetric, and the reflected acoustic wave should have the same shape and size, only with opposite velocity. 17,18 Hence, the pulse at t = 0.5 ought to recover its initial position. Computations consider the flow to be inviscid and are based on several uniform grids with the increasing number of grid cells of N c = 60, 120, 240, 480, 960, and 1920, respectively. A very tiny timestep of Δt = 1 × 10 −5 is used to minimize the temporal integration errors. For the penalization method, we set = 0.01 as suggested in Reference 1 and we found the above timestep is sufficiently small for a simple explicit Euler scheme to solve the penalized density equation in Equation (19). Figure 3 shows the computed profiles of the pressure perturbation at t = 0.5 as the acoustic wave should recover its initial position, centered at x = 0.25. It is shown that for three varying grids, the present variable-extended IBM exhibits perfectly agreeable results with the reference profile, not only in preserving the perturbation amplitude of but also in reducing the phase shift from pulse center x = 0.25. The penalization method, however, gives less accurate results, especially when the grid resolution is coarse below N c = 240. In the penalization method, the spikes inside the solid domain (x > 0.5) are apparent, which is also the case in figure 3(B) of Reference 17. Both methods are able to produce converging solutions as the grid resolution improves. To evaluate the accuracy of both methods, we calculate the errors by comparing the computed wave profile with its initial counterpart within the fluid domain, namely, x ∈ [0, 0.5]. In Figure 4, L 1 and L ∞ error norms of pressure and velocity, respectively, are plotted. As it is shown, although all the error norms follow the first-order scaling with the increasing grid resolution, the present method produces obviously smaller absolute errors than the penalization method. Note that, for both methods, these errors are mostly attributed to how the reflective wall condition at the fluid-solid interface is treated. refined according to the inflow condition and is represented by blocks as shown in Figure 5, including the 61-point circular cylinder zoomed in aside. Note that each block contains a fixed array of 16 × 16 computational cells. Hence, the equivalent spatial resolution would be 1536 × 1024 over the entire domain, except the cylinder boundary with one level more refinement. The cylinder boundary is considered as an adiabatic and nonslip wall so that velocity inside the cylinder is treated by Equation (24). Other boundary conditions consist of the inflow condition at the left side and the outflow condition at the remaining three sides. A fixed CFL number of 0.8 is used. The entire computational domain at t = 0 is initialized using the freestream state except that the velocities inside the solid cylinder are zeros.

C d C l(r.m.s) St
Experimental results 33  In Figure 5 with the vorticity field at the bottom, we can clearly observe the nonsymmetric vortex shedding in the laminar wake behind the cylinder. Time evolutionary lift and drag coefficients C l and C d are plotted in Figure 6, indicating the unsteady vortex shedding tends to be regularly periodic at later stage. In Table 1, we quantitatively compare the present results, which are conveniently sampled and integrated using local quantities at the Lagrangian points of the cylinder, with other numerical and experimental results available. A very satisfactory agreement with reference results in the time-averaged drag coefficient C d and the Strouhal number St can be seen, which implies the capacity of capturing transient behaviors in the flowfield by the present flow solver with the proposed variable-extended IBM. The amplitude of the lift coefficient C l in Figure 6 and the resulted C l(r.m.s) are greater than the reference numerical results which also exhibit a considerably large variation of C l(r.m.s) between themselves. The reason might be attributed to the limited number of sampling points and the interpolation scheme from its surrounding grid cells to each sampling point in statistics.

Case 3: Reactive gas mixture shock tube
We begin to simulate gaseous reactive flows by considering a 1D reactive shock tube to evaluate the performance of the present method compared with the penalization method and traditional "body-fitted" exact reflective wall boundary condition. A reactive gas mixture is placed in a closed tube, in which a shock hits the left-side solid wall boundary and reflects off. A reaction wave occurs at the boundary after a delay of induction, then catches up and merges with the right-moving shock wave to develop a detonation wave. The reactive mixture is characterized by a 2/1/7 molar ratio of with the left and right states being separated in the middle of a 12 cm-long domain. In previous studies, exact wall boundary condition is often applied for the left end of tube using where subscript "g" and "i" denote the ghost cells and their image internal cells, respectively. Instead, in the present IBM and the penalization method, a solid body is immersed within 0 ≤ x ≤ 3 cm to function as the wall and the entire shock tube is translated by 3 cm rightwards. Outflow condition for the right end of tube uses the first-order extrapolation. A 9-species 23-reaction hydrogen/air combustion mechanism 38 is considered here. Results are analyzed at t = 210 μs with N c = 200, 400, and 800 in the fluid domain, respectively, and a fixed CFL number of 0.75. The solid domain, if involved, uses equal grid spacings with the fluid domain. In Figure 7, it can be seen that using the penalized nonslip velocity condition results in slower propagation of the detonation wave at all the three grids. We have attributed such inaccuracy to the loss of total energy at the wall by the penalization method, as mentioned in previous section. By contrast, the present variable-extended method yields much agreeable detonation waves in terms of the location of the reacting front, the velocity peak value and the overall profile of the reflected waves, against the exact wall boundary condition solution. Grid convergence of the present method is especially shown to validate the numerical accuracy of coping with unsteady reactive flows.

Case 4: Shock-induced detonation of Lehr's experiments
We continue to simulate the gaseous reactive flow from Lehr's benchmark experiments, 39 by considering a spherical projectile with radius r = 7.5 mm flying past the premixed reactive gas mixture of H 2 /O 2 /N 2 in three scenarios: (1) super-detonation case. 40,41 The surrounding reactive gas mixture is featured by the molar ratio of H 2 :O 2 :N 2 being 2:1:3.76, static pressure of p ∞ = 320 mmHg and static temperature of T ∞ = 293 K. The projectile velocity is u ∞ = 2605 m/s, corresponding to a supersonic flight Mach number of M = 6.46. (2) subdetonation case. 40 The surrounding reactive gas mixture is featured by the molar ratio of H 2 :O 2 :N 2 being 2:1:0, static pressure of p ∞ = 186 mmHg and static temperature of T ∞ = 293 K. The projectile velocity is u ∞ = 1892 m/s, corresponding to M = 3.55. (3) trans-detonation case. [41][42][43][44] The surrounding reactive gas mixture is the same as that of the superdetonation case. The projectile velocity is lower, that is, u ∞ = 1804 m/s, corresponding to M = 4.48.
At t = 0, the computational grid is adaptively refined according to the initial condition and is represented by blocks as shown in Figure 8. Freestream inflow condition is set for the left boundary of the fluid domain. Boundary conditions are easily imposed for supersonic flows by using the first-order extrapolation. For the bottom side, axisymmetric boundary condition should be applied and additional source terms accounting for the 2D axisymmetric model are given by The spherical projectile surface is considered as an adiabatic wall so that Neumann boundary condition for temperature is applied. As the effect of transport properties is known to be negligible, 42 concerned computations are performed on a basis of the inviscid assumption. A fixed CFL number of 0.75 is used.

F I G U R E 8 (Case 4) Computational domain including the
flow domain and the solid spherical projectile, with the adaptively refined grid of blocks at t = 0 (left) (note that for the superdetonation case, each block contains 16 × 16 grid cells), and the computed steady superdetonation temperature field at t = 4 × 10 −5 s (right), upon which experimental measurements of the shock and combustion front are marked by red circle symbols and black square symbols, respectively. The Evans' hydrogen/air combustion mechanism (seven species/eight reactions) is used 45 Considering the superdetonation first, in Figure 8, the right figure depicts the fully developed steady-state temperature distribution around the spherical head. We can see that near the stagnation streamline, the reacting front of the superdetonation merges with the detached bow shock, because the flow past the shock decelerates so quickly that the high temperature and pressure environment would lead to immediate ignition of the reactive mixture. In the circular direction upwards, bifurcation between the reacting front and the bow shock arises since the induction time/length of ignition tends to be longer. Experimental measurement of the two discontinuity positions is marked upon the present numerical result, matching excellently with each other.
The subdetonation case is described by the computed steady temperature distribution in Figure 9. As it is shown, the reacting zone close to the spherical head is much thinner than that of the superdetonation case, being apart from the detached shock. In comparison with the experimental shadowgraph, it is readily to see that the present numerical result including the shock and reacting front positions and curves, no matter ahead of the projectile or away from the projectile, agrees well with its reference.
The transdetonation case involves an interesting phenomenon of the oscillating reacting front. As shown in Figure 10 (top), the instantaneous temperature fields exhibit obvious unsteady combustion with a noncircular reacting front, which belongs to the regular regime 44 in that the flow oscillates with high frequency and low amplitude. Of all the Lehr's experiments, the Mach 4.48 case gives a measured frequency of 425 kHz. In the present computation, by analyzing the long-time history of pressure at the stagnation point in Figure 10 (bottom left), the calculated frequency of regular cycles by DFT analysis is about 395 kHz for the computation with a block size of 32 × 4 cells. By increasing the y-direction spatial resolution with a block size of 32 × 8 cells, the calculated frequency increases to 410 kHz. Other studies 43,44 also suggest an improved numerical accuracy of the oscillation frequency by increasing the grid resolution. Back to Figure 10 (top), it is also interesting to observe that 1) a small strip of isolated combustion area forms in front of the main detonation wave at t = 57.79 μs, 2) it then develops upwards and merges with the main combustion area at t = 58.75 μs and 3) the merged reacting front continues to move along the circular direction at t = 59.71 μs. In this process, we can also find that the thickness of the main combustion area along the stagnation streamline increases at the beginning to reach a peak value and then decreases. These three stages, from Figure 10 (bottom right), in fact correspond to the minimum, approximately the maximum and the mean value, respectively, of the oscillating pressure at the stagnation point in a cycle. As a result, the pressure oscillation at the stagnation point is believed to be highly related to the repeated stretch and shrink of the combustion area attached to the projectile head, which challenges the accuracy and robustness of the numerical method in use and verifies the present variable -extended IBM.
F I G U R E 9 (Case 4) Comparison of the computed steady subdetonation temperature field (mirrored against the x-axis) with the experimental shadowgraph. 39,40 Note that for the subdetonation case, each block contains 16 × 4 grid cells to reduce computational costs. The Jachiomowski's hydrogen/air combustion mechanism (9 species/19 reactions) is used 42 F I G U R E 10 (Case 4) The unsteady transdetonation temperature fields at three time instants (top), the histories of pressure at the stagnation point over time using two grids of blocks (bottom left) and the zoomed-in pressure oscillation cycle (bottom right). The Jachiomowski's mechanism (9 species/19 reactions) is also used here

CONCLUSIONS
In this article, a variable-extended IBM has been proposed to simulate compressible gaseous reactive flows past solid bodies. In order to alleviate the numerical stiffness of the penalized continuity equation as well as to facilitate the implementation of the Neumann boundary condition in BPM, an extending step from the fluid to solid domain within a narrow band is developed for relevant primitive variables including but not limited to velocity, pressure and temperature. With the help of an EoS, density of ghost cells inside the solid domain can be determined using the extended steady-state pressure and temperature, instead of solving the stiff penalized continuity equation directly. The adiabatic nonslip or slip wall condition can be thereby easily imposed. Meanwhile, the advantage of the BPM in avoiding the local geometric complexity of the solid-fluid interface is maintained, without considering the cut cell types, intersection points and image points. Combined with an adaptive Cartesian grid-based flow solver for multispecies compressible Navier-Stokes equations with finite-rate nonequilibrium chemistry, the present variable-extended IBM has been tested by 1D and 2D numerical examples, and exhibits satisfactory performance in capturing (1) the reflective nature of fluid-solid interface and (2) the nonreactive or reactive, weakly compressible or high Mach number, steady or unsteady flowfields in a simple and robust manner. Besides, treating the velocity carefully in the ghost cells inside the solid domain helps to conserve the total energy and results in more accurate propagation of the detonation front in energy-sensitive reactive flows. Some future work concerning the solid body dynamics and more complicated multiphase flows can be conducted.

ACKNOWLEDGMENT
The financial support from the EU Marie Skłodowska-Curie Innovative Training Networks (ITN-ETN) (Project ID: 675528-IPPAD-H2020-MSCA-ITN-2015) for the first author is gratefully acknowledged. Open access funding enabled and organized by Projekt DEAL.

DATA AVAILABILITY STATEMENT
Some or all data, models, or code generated or used in this study are available from the corresponding author by request.