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 body-conformal 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 non-conforming 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 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. According to the work of Cui et al. [11], the IBM can be classified into two categories, namely diffused-interface 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 [10,12] and a source term [13] to the surrounding Cartesian grid points via a delta function [14]. After that, these two terms are added to the momentum equation and energy equation, respectively. Our previous work [15] developed an implicit diffused-interface IBM based on variable correction for simulating compressible viscous flows over stationary and moving bodies. 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 [12,13].
However, Sotiropoulos and Yang [16] 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 regard as a discrete forcing approach, in which the boundary is precisely tracked. It includes cut-cell method [17,18], Cartesian IB method [11,19], 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, 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. [20] 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,20]. 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. Shuvayan et al. [19] 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,11] to reconstruct the values at IB points. Khalili et al. [9] used a bilinear interpolation for two-dimensional problems, and Mittal et al. [20] 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 rhoCentralFoam 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 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, t e is the total energy that includes kinetic energy and internal energy: where e is the internal energy, v C 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 is defined as
where U ∞ is the freestream velocity, 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 B p 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 ( c U , c X ) 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.

( )
,t Γ n = n X 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 cell-centered finite-volume method. High-speed compressible flows can be simulated by the solver rhoCentralFoam, which is a semi-implicit segregated density-based solver with the Kurganov-Tadmor divergence scheme [21]. More details about the algorithmic structure of rhoCentralFoam solver can be found in Ref. [13]. In our recent work [15], 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 no-penetration boundary condition should be used. The 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 ( ) c t U 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 / 0 p n ∂ ∂ = can be applied on the boundary. But it should be reformulated for a moving body. Following the work of Udaykumar et al. [22], 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. / D Dt u can be approximated directly from the known boundary velocities and this obviates the approximation of the convective term [8,22].
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 [15], 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 where ( ) The unknown coefficients C can be computed from the four surrounding points with their flow information: If i P 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: 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. Fig. 2  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 x y x sign n x y x y x y sign n y x y x sign n x y sign n y x y represents the coordinate of the i P , sign is the symbolic function, and , x y ∆ ∆ are the mesh steps. If i P 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

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 BI p can be achieved by using Eq. (10). Then Eq. (7) can be rewritten as: where i s ∆ 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 time level n+1 is given as follows: 1) For a passively moving body, solve its motion and get c U and c X by using Eqs. (5) and (6); Otherwise, c U and c X are given in advance.  In current simulations, the rectangular computational domain is ( )  / p p , temperature ratio 2 1 / T T , 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.

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 (1) (3) (4) and 4 (2). For the case of 0.85, 1 Ma α = =  , however, it is indicated that such problem is not easy to be handled [24,25]. The strong shock appears on both the upper surface and the lower surface, as shown in Figs. 4(3) and 4(4). (1) (2)  Based on two problems above, the sharp-interface IBM in this study has been proved to be able to deal with stationary boundary problems accurately. Then the piston is given a sudden velocity with 2.0 Ma = . The same problem has also been simulated in previous studies [17,18,27], and it has an analytical solution. 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.

Piston moving with supersonic velocity
Using the finest mesh, Fig. 8  (1)  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 1 0.2 m m × , as illustrated in Fig. 9. At the initial moment, a shock wave with Mach 3 is located at 0.08 , which is traveling from left to right. The pre-shock state, which is on the right of the shock, is at rest with 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 [28,29]. , respectively. Obviously, with the refinement of the grid, the numerical result will gradually tend to an accurate value, both in Shyue's work [29] and present work. Using the same mesh step, the present results are in good agreement with those of Shyue [29]. As shown in Fig. 10 .
Similarly, Fig. 11 shows the current results at the same instants. From two figures, it is found 26 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. (1) (2) Fig. 10 Pressure contours computed by Shyue [29] at different time: (1) 0.1641 t s = , and (2) 0.30085 t s = . (1) (2)

Conclusions
A robust sharp-interface immersed boundary method is proposed in this work for 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.