A sharp-interface immersed boundary method for simulating high-speed compressible inviscid flows

This paper presents a robust sharp-interface immersed boundary method for simulating inviscid compressible flows over stationary and moving bodies. The flow field is governed by Euler equations, which are solved by using the open source library OpenFOAM. Discontinuities such as those introduced by shock waves are captured by using Kurganov and Tadmor divergence scheme. Wall-slip boundary conditions are enforced at the boundary of body through reconstructing flow variables at some ghost points. Their values are obtained indirectly by interpolating from their mirror points. A bilinear interpolation is employed to determine the variables at the mirror points from boundary conditions and flow conditions around the boundary. To validate the efficiency and accuracy of this method for simulation of high-speed inviscid compressible flows, four cases have been simulated as follows: supersonic flow over a 15 ° angle wedge, transonic flow past a stationary airfoil, a piston moving with supersonic velocity in a shock tube and a rigid circular cylinder lift-off from a flat surface triggered by a shock wave. Compared to the exact analytical solutions or the results in literature, good agreement can be achieved.


Introduction
Inviscid flows are usually used as flow field conditions for fast transient problems, such as the impact of shock waves onto a structure [1][2][3][4] and the internal flow field arising within a solid rocket motor [5]. To treat such problems, the conventional body-conformal grid methods [6,7] can be utilized. Although it is direct to enforce boundary conditions on the fluid-structure interface with a bodyconformal grid method, the complex-shaped structure increases too much work on the grid generation process. Moreover, a re-meshing algorithm is required to adapt a moving solid boundary which consumes more computing time and increases the algorithmic complexity of these body-conformal grid methods. For using a nonconforming Cartesian grid, on the other hand, the immersed boundary method (IBM) has gained attention as an appropriate method for tackling flows with complex-shaped stationary or moving boundaries during the last two decades. Applying the advantage of Cartesian grid, in particular, some high-order algorithms for complicated moving problems have been developed [8,9].
It is well known that the IBM was firstly proposed by Peskin [10] to simulate the blood flow in the heart. Following the pioneer work of Peskin, numerous scientists have contributed to improve the accuracy and efficiency of the IBM. Huang et al. [11] introduced the variety of IBM fundamentals and assessed the latest progresses especially the strategies to address the challenges and the applications of the IBM. According to the work of Cui et al. [12], the IBM can be classified into two categories, namely diffusedinterface methods and sharp-interface methods based on the representation of the fluid-structure interface. The diffused-interface IBM can be regarded as a continuous forcing approach. The boundary is smeared by distributing a forcing term [11,13] and a source term [14] to the surrounding Cartesian grid points via a delta function [15]. After that, these two terms are added to the momentum equation and energy equation, respectively. Our previous work [16] developed an implicit diffused-interface IBM based on variable correction for simulating compressible viscous flows over stationary and moving bodies. Furthermore, Wang et al. [17] applied a diffused-interfaced IBM to simulate compressible multiphase flows. The significant advantage of diffused-interface IBM is that it is formulated independent of the spatial discretization, and it simplifies the implementation into an existing fluid solver, just like [13,14]. However, Sotiropoulos and Yang [18] indicated that the classic diffused-interface IBM is difficult to simulate fluid-structure interaction problems, because the inherent stiffness of the forcing term may induce numerical instabilities and spurious oscillations.
On the other hand, the sharp-interface IBM can be regarded as a discrete forcing approach, in which the boundary is precisely tracked. It includes cut-cell method [19,20], Cartesian IB method [12,21], and ghost-cell method [8,9]. Among these methods, the cut-cell method can obtain the clearest interface. The reason is that the boundary divides grid cells into two sub-cells for the solid phase and fluid phase, respectively. Therefore, the conservation of mass, momentum and energy can be strictly enforced into the fluid phase. However, the complicated cell reshaping procedure causes difficulties in simulating complex moving body problems. In the Cartesian IB method, the fluid points with at least one neighboring point inside the solid body are called IB points or forcing points. Reconstructing the flow variables at these points and enforcing the boundary conditions can ensure the accuracy of the flow field. The ghost-cell method is similar to the Cartesian IB method. The solid points with at least one neighboring point inside the fluid domain are called ghost points. The reconstruction procedure is done at these ghost points. Tran and Plourde [5] used both the Cartesian IB method and the ghost-cell method to solve internal flows. Although the Cartesian IB method and the ghost-cell method are less accurate than the cut-cell method because of their implicit representation of the solid boundary, the point recognition procedure is easier than the cut-cell procedure, and the flux calculation around the immersed boundary is also not necessary. Mittal et al. [22] have shown the large potential of ghost-cell method for simulation of highly complex moving body problems.
In the ghost-cell method, the flow variables at the ghost points cannot be interpolated from the fluid field or the boundary conditions directly. Thus, an image point can be considered, which is the mirror of ghost point along the normal direction to the boundary and is always in the fluid domain. A body intercept point is the midway point between ghost point and image point, which must be on the boundary. Therefore, boundary conditions can be enforced on body intercept points directly. The flow variables at ghost points are interpolated from corresponding image points and body intercept points [5,8,9,22]. A significant issue of ghost-cell method is the accuracy of the interpolation procedure at image points from corresponding body intercept points and fluid points. Brahmachary et al. [21] used an inverse distance weighting interpolation method to reconstruct the flow variables. This scheme preserves the local maxima and minima, which leads to a smooth reconstruction and is more stable than a polynomial reconstruction. A linear interpolation method was used in [5,12] to reconstruct the values at IB points. Khalili et al. [9] used a bilinear interpolation for two-dimensional problems, and Mittal et al. [22] used a trilinear interpolation for three-dimensional problems similarly. Higher-order polynomials are required for more accurate results, but they often lead to numerical instabilities, and more information is needed to determine the polynomials. Qu et al. [8] listed four polynomials for the interpolation: a linear, an incomplete quadratic (bilinear), a complete quadratic and a cubic. They found that higher-order polynomials (e.g., a cubic polynomial) required more interpolation points and the computational time may significantly increase. Furthermore, they revealed that the complete quadratic polynomial does not significantly improve the accuracy of their method but requires more interpolation points than the incomplete quadratic polynomial. It is shown that bilinear interpolation may be an expected scheme to reconstruct the flow variables.
In this work, we present a robust IBM to simulate high-speed compressible inviscid flows over stationary and moving bodies. The Euler equations are discretized on a Cartesian grid and solved by a segregated density-based solver, namely rho-CentralFoam with Kurganov and Tadmor divergence scheme in the open source library OpenFOAM. The sharp-interface IBM is adopted, which is a significant issue in this work, because the exact location of shock waves has to be predicted. Khalili et al. [9] developed a ghost-cell IBM to simulate viscous compressible flows, and a bilinear interpolation was utilized to reconstruct the flow variables. However, it is noted that only low Mach number flows were considered in their study and the obtained conclusions were similar to that of incompressible flows. In this study, the similar ghost-cell IBM is adopted for enforcing the boundary conditions in simulating high-speed compressible inviscid flows. For a moving rigid object, moreover, the Newton's second law of motion is utilized to solve strong coupling between fluid and structural within each time step. Because of inviscid flows, there is no viscous stress or torque on solid surface and the rotational motion is not involved. To validate the proposed method combined with OpenFOAM, two stationary cases are simulated: supersonic flow over a 15°angle wedge and transonic flow past a stationary airfoil. Then a moving piston and a lifting-off circular cylinder are simulated to validate the current method for handling moving body problems. All the results are compared well with the exact analytical solutions or the results in literature.
This paper is organized as follows. In Section 2, governing equations for fluid flows and motions of rigid bodies are described. In Section 3, numerical methods are presented, i.e., the governing equations in OpenFOAM and the ghost-cell IBM. In Section 4, four cases are simulated to validate the efficiency and accuracy of this method. In Section 5, conclusions are summarized.

Governing equations
For flows over a stationary body, only the flow field should be resolved. But for flows over a moving body, both the flow field and the motion of body should be resolved. Moreover, the fluid is strongly affected by the existence of body, and the body motions are coupled with the fluid flows.

Governing equations for fluid flows
For inviscid compressible flows, the Euler equations are solved, which can be written as follows: where ρ is the fluid density, u is the fluid velocity, and p is the pressure. In addition, e t is the total energy that includes kinetic energy and internal energy: where e is the internal energy, C v is the heat capacity at constant volume, and T is the fluid temperature. As the perfect gas is considered here, the equation of state p = ρRT is adopted, where R is the specific gas constant.
The Mach number Ma is defined as ocity, T ∞ is the freestream temperature, and γ is the specific heat ratio, which is kept constant at 1.4 in present work. The pressure coefficient is defined as where p B is the pressure on the boundary of body, p ∞ is the freestream pressure, and ρ ∞ is the freestream density.

Governing equations for a passively moving rigid body
For a moving rigid body, only the velocity and the position vector of the mass center (U c , X c ) are considered in this work. Due to the inviscid flows, the rotational motion of the body is not involved since there is no existence of viscous stress or torque. According to the Newton's second law, the translational motion of a passively moving rigid body can be solved by the following equations: where m is the mass of the body, and F f is the resultant force imparted by the fluid to the body, which can be given by: where f f represents the force per unit area on the boundary of body by the fluid, which only includes the pressure interpolated from the flow fields, and Γ represents the boundary of body. n = n(X Γ , t) is the unit normal vector from the boundary into the flow field, where X Γ represents the position of the point on the boundary.

Basic solver in OpenFOAM
As mentioned in Introduction, an open source code, i.e. OpenFOAM, is used to solve the governing equations of flow field. The codes in OpenFOAM are based on a cell-centered finitevolume method. High-speed compressible flows can be simulated by the solver rhoCentral-Foam, which is a semi-implicit segregated density-based solver with the Kurganov-Tadmor divergence scheme [23]. More details about the algorithmic structure of rhoCentralFoam solver can be found in Ref. [14]. In our recent work [16], this solver has been combined with a diffused-interface IBM to simulate viscous compressible flows. To extend this solver to deal with inviscid compressible flows, only the viscous terms of the governing equations are ignored.

Boundary conditions
Classical boundary conditions including subsonic or supersonic inlet, subsonic or supersonic outlet and walls are applied at the external boundaries of the computational domain. For the inviscid flows, no boundary layer should develop in the vicinity of the boundary of body. When a stationary body is considered, the traditional slip boundary condition can be used on the boundary of body [5]. But when a moving body is considered, the nopenetration boundary condition should be used. The velocity vector u(X Γ , t) on the boundary Γ is decomposed into the normal and tangential velocity components: where τ(X Γ , t) is the unit tangential vector at X Γ , which is perpendicular to the unit normal vector. Considering the no-penetration boundary condition, these components are then written as: where the velocity of the body mass center U c (t) has been introduced in the previous section and it is independent of X Γ for a rigid body. For a stationary body, the pressure boundary condition ∂p/∂n = 0 can be applied on the boundary. But it should be reformulated for a moving body. Following the work of Udaykumar et al. [24], the pressure gradient normal to the boundary is obtained by projecting the momentum equation into the boundary normal direction, which results in: where D/Dt is the material derivative. Du/Dt can be approximated directly from the known boundary velocities and this obviates the approximation of the convective term [8,24]. In addition, an adiabatic boundary condition is applied at the boundary. The temperature gradient normal to the surface is written as: The boundary conditions above can be classified into two categories, Dirichlet boundary condition (Eq. (10)) and Neumann boundary conditions (Eqs. (11), (12), (13)). In our previous work [16], a variable correction-based IBM is developed to correct the variables with the Dirichlet boundary condition, including the no-slip boundary condition and the isothermal boundary condition. But it is difficult to handle the Neumann boundary conditions. In this work, both two boundary conditions need to be enforced into the proposed method.

Sharp-interface immersed boundary method
As one of sharp-interface immersed boundary methods, the ghost-cell IBM is employed to impose the boundary conditions on the boundary of body in this work. According to the location of boundary, the computational grid nodes lying inside the body are identified as the solid points, and the nodes lying outside the body are identified as the fluid points. Considering the rhoCentralFoam solver employed, one layer of solid points adjacent to the immersed boundary are called ghost points (GPs), as illustrated in Fig. 1. The interpolation procedure cannot be conducted on these GPs directly. Thus, an image point (IP) is required, which is the mirror of GP along the normal direction to the boundary. Thus, a boundary intercept (BI) point is defined as the midway point between GP and IP. The unit normal vector n is also marked on the BI point in Fig. 1.
After identifying different kinds of points, the interpolation procedure can be conducted on the IP. In Ref. [5], a liner interpolation was used and it followed the formulation ϕ = c 1 x + c 2 y + c 3 . Three coefficients (c i , i = 1, 2, 3) should be computed by using three points with known values. The IP should be surrounded by these three points to maintain good stability. However, each cell in the Cartesian grid is rectangular, and the Fig. 1 Different kinds of points on a Cartesian grid IP must be in a rectangular domain, as illustrated in Fig. 1. The four vertices of the rectangular domain can be exactly regarded as known points. The liner interpolation should choose three of them to surround the IP, which may increase the complexity of the algorithm. Alternatively, a bilinear interpolation following the formulation ϕ = c 0 + c 1 x + c 2 y + c 3 xy [8] is adopted in this work, and four unknown coefficients (c i , i = 0, 1, 2, 3) are determined by the four surrounding points. The accuracy of this interpolation scheme is second order, which is accurate enough for the solver in OpenFOAM. Thus, this interpolation method is simpler but more accurate than the linear interpolation algorithm. A generic variable ϕ can be expressed as: where C = [c 0 , c 1 , c 2 , c 3 ] T is a coefficient vector, and ε ¼ 1; x; y; xy ½ T is a coordinate vector that only depends on a local coordinate system (x, y). As mentioned above, the coefficient vector C should be determined by four surrounding points. However, if the surrounding point is just one GP, its corresponding BI would be used. Then the interpolation domain changes from a rectangle to the blue area in Fig. 1. The four surrounding points are marked as P i , i = 0, 1, 2, 3, which can be a fluid point or a BI point. Due to the particularity of the BI point, Khalili et al. [9] considered multiple situations for interpolation according to the number of the BI points in the interpolation domain. In the current method, however, such complicated procedure is discarded, and only a judgment step is required to determine the information type of a surrounding point. The flow variables ϕ on the fluid point can be obtained from the flow field, while the values on the BI point can be obtained by using the boundary conditions. The Dirichlet-type boundary condition (Eq. (10)) is enforced on the BI point just like the fluid point. But for the Neumann-type boundary conditions (Eqs. (11), (12), (13)), some modifications are required: where (x, y) is the BI point coordinate, and n = (n x , n y ) is the unit normal vector at the BI point.
The unknown coefficients C can be computed from the four surrounding points with their flow information: where A ¼ ε 0 ; ε 1 ; ε 2 ; ε 3 ½ T is a 4 × 4 matrix depending on the coordinate information of the four surrounding points, and Φ ¼ ½ϕ 0 ; ϕ 1 ; ϕ 2 ; ϕ 3 T is a generic variable vector. ε i and ϕ i (i = 0, 1, 2, 3) can be obtained by a judgment step, based on the information type of P i : If P i is a BI point and the given boundary condition is the Neumann-type boundary condition, Eq. (15) will be used.
Then the flow variables at the image point ϕ IP can be obtained as: After obtaining the values at the IP, the flow variables at the GP can be determined by using the boundary conditions on the immersed boundary. A linear interpolation along the normal is employed. Different types of boundary conditions need different formulas to calculate the values of flow variables at the GP: where Δl is the length of the normal probe from the GP to the IP.

Fresh points for a moving body
For the case of stationary body, the GPs can be identified initially. But for the case of moving body, the GPs have to be updated at each time step. Figure 2 shows the process of immersed boundary moving from time level n to time level n + 1. As a result, a fluid point at time level n may become a newly emerged GP at time level n + 1 when the immersed boundary intrudes into the fluid domain, or a GP at time level n may become a newly emerged fluid point (called "fresh point", i.e., FP) at time level n + 1 when the immersed boundary withdraws from the fluid domain. For the first situation, the newly emerged GP can be regarded as a normal GP and the reconstruction procedure is as usual. For the second situation, however, the FP is not possible to devise a straightforward temporal discretization of the governing equations, since it lacks a settled time history. Tran and Plourde [5] reconstructed the flow variables both on ghost point and forcing point. The forcing points are the points with at least one neighboring point inside the body. Since the FPs are included in the forcing points, this method seems to be workable for the moving boundary problems. Nevertheless, it is found that the FPs do not always appear, and this method also has to consume too much computing time.
In the present method, the reconstruction procedure on the FP is similar to that of GP. The interpolation domain changes from a rectangle to the blue area in Fig. 2, as the FP is shifted to its corresponding BI point. The other three surrounding points are marked as P i , i = 1, 2, 3, which can be confirmed by: where (x i , y i ) represents the coordinate of the P i , sign is the symbolic function, and Δx, Δy are the mesh steps. If P i is a FP, it will be also shifted to its corresponding BI point.
And this manipulation is similar to that of GP. The flow variables are directly interpolated to the FP: At each time step, the layer of FP is required to be one layer depth at least to guarantee the success of interpolation step. Considering the motion of immersed boundary and its velocity, the time step should be constrained as Δt < min(Δx, Δy)/|U c |.
It should be clarified that the current interpolation procedure has some general similarities with the method proposed by Mittal et al. [22]. However, their interpolation stencil was always in a rectangular domain and the surrounding points would contain other GPs. As a result, a point Gauss-Seidel method was used for these coupled GPs.
In the current method, if the surrounding point is just one GP, its corresponding BI point will replace it. Consequently, the surrounding points only include fluid points or BI points in the current interpolation stencil, as illustrated in Fig. 1. Then, the flow variables at the GP can be computed easier. Another advantage of the current method is that the BI points are included in the surrounding points and the boundary information is enforced into the interpolation procedure. The second difference is the reconstruction procedure of FPs. Mittal et al. [22] used the same method as GPs to calculate FPs. When the FP is in the fluid domain, the current method directly uses surrounding points to perform the interpolation to obtain the flow variables at the FP. Therefore, the calculation process is more concise.

Computing the forces on the body
For inviscid flows, only the pressure acting on the surface of body is required to be computed. The BI point is treated as an interpolation point, and the pressure p BI can be achieved by using Eq. (10). Then Eq. (7) can be rewritten as: where Δs i is the arc length of the immersed boundary around the BI point. The pressure coefficient at the BI point can be expressed as:

Solution procedure
To summary, the solution procedure of the current algorithm from time level n to n + 1 is given as follows: 1) For a passively moving body, solve its motion and get U c , X c by using Eqs. (5) and (6); Otherwise, U c and X c are given in advance. the density ρ at GPs is updated by using the equation of state; 5) Solve the governing equations of flow field, and compute the flow variables by using OpenFOAM; 6) Calculate the resultant force F f by using Eq. (24); 7) Repeat Steps 1 to 6 until convergence is reached, and obtain ρ n + 1 , u n + 1 , p n + 1 and T n + 1 .

Supersonic oblique shock
To test the accuracy of the present method, a supersonic flow over a 15°angle wedge is simulated. This flow can produce a supersonic oblique shock, which has an exact analytical solution. Two Mach numbers are considered, i.e., Ma = 3.0 and Ma = 5.0. Supersonic inlet and zero-gradient outlet are adopted. Initial conditions are chosen with p = 1Pa and T = 1K. In current simulations, the rectangular computational domain is (1.5m, 1.0m). A very fine uniform Cartesian grid is used with the mesh size of 1500 × 1000, and then the mesh spacing is Δh = 0.001m. At Ma = 3.0, the pressure contours are exhibited in Fig. 3. The slope starts at x = 0.5m, where the supersonic flow forms an oblique shock wave. The wave surface of the oblique shock is very clear, which is attributed to the good ability of the current method to handle the supersonic flow over a body. Table 1 lists the theoretical and numerical results from pre-shock state (subscript 1) and post-shock state (subscript 2): post-shock Mach number Ma 2 , pressure ratio p 2 /p 1 , temperature ratio T 2 /T 1 , and shock-angle that is the angle between the oblique shock and the horizontal plane. The present results are in agreement with the exact solution, and the accuracy of the current method is well displayed.
To assess the convergence rate of the current method, the case at Ma = 3.0 is simulated with different meshes. Four uniform grids of 75 × 50, 150 × 100, 300 × 200 and 1500 × 1000 are used. The corresponding grid steps are Δh = 0.02m, Δh = 0.01m, Δh = 0.005m and Δh = 0.001m, respectively. As illustrated in Table 1, the exact solution of shock-angle is θ = 32.24 ∘ . The error of shock-angle e θ versus grid steps in the log scale is plotted in Fig. 4. From the figure, it is found that the slope of the line is about 1.7. Although the accuracy of the interpolation scheme is second order in the current IBM, some discrete schemes in the rhoCentralFoam solver are only first-order. Therefore, the overall accuracy of the current solver is less than the second-order.

Transonic flow around a stationary airfoil
To verify the robustness of the present method, a flow around an irregular-shaped object is simulated in this section. Here, a NACA0012 airfoil is considered. Two transonic flows are selected, in which the free-stream Mach number and angle of attack are: (a) Ma = 0.8, α = 1.25 ∘ and (b) Ma = 0.85, α = 1 ∘ . The same problem has been simulated in previous work [25][26][27]. The computational domain is 28c × 20c with the mesh size of 511 × 251, where c is the chord length of airfoil, and the airfoil is located at (10c, 10c). The region around the airfoil is 1.2c × 0.2c, and a fine uniform mesh with the spacing of h = 0.005c is adopted.
The pressure and Mach number contours are presented in Fig. 5. For the case of Ma = 0.8, α = 1.25 ∘ , there are a relatively strong shock on the upper surface of the airfoil and a weak shock on the lower surface of the airfoil, as illustrated in Fig. 5(1) and 5(2).  For the case of Ma = 0.85, α = 1 ∘ , however, it is indicated that such a problem is not easy to be handled [25,26]. The strong shock appears on both the upper surface and the lower surface, as shown in Fig. 5(3) and 5(4). Figure 6 shows the pressure coefficient −C p on the upper and lower airfoil surfaces. From the figure, it can be found that the present results compare well with the results of previous work [27,28]. At Ma = 0.8, α = 1.25 ∘ , the shock positions on the upper and lower airfoil surfaces are accurately predicted. Although the case (b) is difficult to simulate, the current method still can accurately predict the position of the upper and lower surface shock waves.
Based on the two problems above, the sharp-interface IBM in this study has been proved to be able to deal with stationary boundary problems accurately.

Piston moving with supersonic velocity
To further validate the current method, the moving body problems will be simulated. The first one is a piston moving at a constant speed of Ma = 2.0 in a quiescent fluid. The movement sketch is shown in Fig. 7. The thickness of the piston is L, and it is located in the middle of a shock tube initially. The shock tube has a length of 128L and a width of 4L. In the whole computational domain, the parameters of the initial quiescent fluid are u = v = 0 m/s, p = 1Pa, and T = 1K. Then the piston is given a sudden velocity with Ma = 2.0. The same problem has also been simulated in previous studies [19,20,29], and it has an analytical solution.
The computational domain with the size of 128L × 4L is discretized by using three different uniform grids with 512 × 16, 1280 × 40 and 2560 × 80. These three grids   Figure 8 shows the pressure distribution along x-direction with different mesh steps. The piston is represented by a blue strip in the figure. p ref = 1Pa is the reference pressure. Obviously, the result from fine mesh is closer to the exact solution than that from coarse mesh, especially at the corners of the curve. However, even for the coarsest mesh, the calculated result is generally in good agreement with the exact solution. Thus, it means that the current method has good robustness, and the satisfactory result can be obtained even at a relatively coarse grid.  Using the finest mesh, Fig. 9 shows the pressure ratio p/p ref , density ratio ρ/ρ ref , and Mach number ratio Ma/Ma ref distribution along x-direction of shock tube. All the results compare well with the analytical solutions. This testing problem verifies that the current method can accurately solve the moving body problem.

Lift-off of a cylinder
The second moving body problem to be simulated is a rigid cylinder whose motion is induced by the fluid. This problem takes place in a shock tube, whose domain size is 1m × 0.2m, as illustrated in Fig. 10. At the initial moment, a shock wave with Mach 3 is located at x = 0.08m, which is traveling from left to right. The pre-shock state, which is on the right of the shock, is at rest with p = 1Pa, T = 1K, and u = v = 0 m/s. The flow variables at the post-shock state can be calculated from the total temperature. The radius of rigid cylinder is 0.05 m, which is at rest initially on the floor of the shock tube. The center of the cylinder is initially at (0.15m, 0.05m), and the density of the cylinder is 10.77kg/m 3 . When the Mach 3 shock wave moves to the cylinder position, the cylinder will be lifted due to the action of the fluid. This cylinder lift-off problem has also been numerically studied in previous work [30,31]. When the shock wave reaches the cylinder, the cylinder can affect the flow field significantly to make the flow field asymmetrical. Due to the strong fluidstructure interaction, this asymmetry causes unequal forces on the upper and lower surfaces of the cylinder, and thus the lift-off phenomenon happens. The instant at t = 0.30085s is set as a stopping time of simulation when the cylinder is close to the upper surface of the shock tube. In Table 2, a convergence study of position of the cylinder mass center is presented. The computational domain is discretized by five uniform grids corresponding to the five grid steps of Δh = 1/500m, Δh = 1/1000m, Δh = 1/1600m, Δh = 1/2000m and Δh = 1/2500m, respectively. Obviously, with the refinement of the grid, the numerical result will gradually tend to an accurate value, both in Shyue's work [31] and present work. Using the same mesh step, the present results are in good agreement with those of Shyue [31]. As shown in Fig. 11, Shyue [31] provided the pressure contours at two different moments, i.e., t = 0.1641s and t = 0.30085s. Similarly, Fig. 12 shows the current results at the same instants. From two figures, it is found that the position of shock waves is captured accurately by the present method. This testing problem completely demonstrates the robustness and accuracy of the present sharp-interface IBM.

Conclusions
A robust sharp-interface immersed boundary method is proposed in this work for simulating high-speed inviscid compressible flows over stationary and moving bodies. The Euler equations are discretized on a Cartesian grid and solved by the rhoCentralFoam solver in the platform of OpenFOAM. The interaction of flow and body boundary is handled by using a ghost-cell IBM. A bilinear interpolation method is employed to robustly interpolate the flow variables at  an image point corresponding to those at a ghost point, as well as to satisfy the boundary conditions. The velocity around the body boundary is decomposed into the normal and tangential components. The normal component satisfies the no-penetration boundary condition, which is the Dirichlet-type boundary condition. The tangential component together with the pressure and the temperature satisfies the Neumann-type boundary conditions at the immersed boundary. Finally, the density is reconstructed via the equation of state at ghost points. Additionally, a newly emerged fresh point may appear when the body moves. The flow variables at fresh points are reconstructed by using a method similar to the reconstruction of ghost points. To validate the robustness and accuracy of the proposed method for high-speed inviscid compressible flows over stationary and moving bodies, four simulations are  performed: supersonic flow over a 15°angle wedge, transonic flow past a stationary airfoil, a piston moving with supersonic velocity, and a rigid circular cylinder lift-off. It is shown that the results from the current method are in good agreement with the data in literature. As the future work, the current sharp-interface IBM can be extended to simulate high-speed viscous compressible flows together with the turbulence model for high Reynold number situation.