An immersedboundary method for compressible viscous flow and its application in gaskinetic BGK scheme
Abstract
An immersedboundary (IB) method is proposed and applied in the gaskinetic BGK scheme to simulate incompressible/compressible viscous flow with stationary/moving boundary. In present method the ghostcell technique is adopted to fulfill the boundary condition on the immersed boundary. A novel idea “local boundary determination” is put forward to identify the ghost cells, each of which may have several different ghostcell constructions corresponding to different boundary segments, thus eliminating the singularity of the ghost cell. Furthermore, the socalled “freshcell” problem when implementing the IB method in movingboundary simulation is resolved by a simple extrapolation in time. The method is firstly applied in the gaskinetic BGK scheme to simulate the TaylorCouette flow, where the secondorder spatial accuracy of the method is validated and the “superconvergence” of the BGK scheme is observed. Then the test cases of supersonic flow around a stationary cylinder, incompressible flow around an oscillating cylinder and compressible flow around a moving airfoil are conducted to verify the capability of the present method in simulating compressible flows and handling the moving boundary.
Keywords: immersedboundary method, moving boundary, compressible flow, gaskinetic scheme
1 Introduction
In recent years, rising attention has been paid to a class of nonboundaryconforming grid methods called the immersedboundary (IB) method for its convenience in handling problems with complex as well as moving boundaries. In the IB method, the grid points don’t have to coincide with the solid boundary and a simple uniform or nonuniform Cartesian grid can be used, which will make the grid generation much easier than the conventional boundaryconforming grid method and guarantee the grid quality for even an extremely complex boundary. Without grid transformation and regeneration, the IB method is very applicable to moving boundary problems with relatively simple solution procedure and significantly lower computational cost compared to the conventional boundaryconforming grid method, especially for the problem involving multibody movement, such as the problem concerns the interaction between fluid and fragments of aircraft, where the laborious grid regeneration or transformation can drive one crazy.
According to Mittal and Iaccarino [1], most of the existing IB methods can be categorized into two groups, i.e. the continuous forcing approach and the direct (or discrete) forcing approach. In the continuous forcing approach, a forcing term is added to the continuous governing equations before they are discretized. The force is determined from the boundary condition on the Lagrangian boundary point and the discrete Dirac delta is usually used to link the Lagrangian and the Eulerian variables. The major advantage for the continuous forcing approach is that it can be easily implemented in an existing numerical scheme to deal with problems with complex or moving boundaries. The main drawback is the low spatial accuracy (generally firstorder local spatial accuracy) due to the use of the discrete delta function. Besides, this class of IB methods is rarely used in the compressible problem.
In the direct forcing approach, the boundary condition is imposed by directly constructing variables on the grid points or cells near the immersed boundary involved in the calculation stencil of the numerical scheme. In contrast to the continuous forcing approach, for this category of IB methods the implementation is more complex and the inclusion of boundary motion is more difficult. But the advantage is that it allows a sharp representation of the immersed boundary and secondorder local spatial accuracy can be obtained for some of the approaches. The pioneer work for this category of methods was done by MohdYusof [2] and Fadlun et al. [3]. De Palma et al. [4], Ghias et al. [5] and de Tullio et al. [6] extended the direct forcing approach to the viscous compressible problem. Among these approaches, quite a part of them are based on the idea of “ghost cell”, such as the methods of Gibou et al. [7], Majumdar et al. [8], Tseng and Ferziger [9], Ghias et al. [5], Mittal et al. [10] and others. One problem for the ghostcell approach is that when dealing with a convex boundary, a ghost cell may be adjacent to more than one fluid cells and there may be several boundary segments lie in more than one directions of the ghost cell. In this singular case, most of the past methods will construct the ghost cell according to the nearest boundary segment. However, it is not rational to assume a unique construction for the ghost cell in such a case. In fact, in our work, some suspicious evidence for this singularity problem can be observed (detailed in Subsection 3.1).
In this paper, we propose an IB method based on the ghostcell approach to simulate the incompressible/compressible viscous flow with stationary/moving boundary. The ghost cell is identified and constructed by an idea of “local boundary determination”, which allows more than one construction for a sole ghost cell, thus no singularity occurs for the ghost cell. The conception of the “image point” presented in some previous methods [8, 9, 5] is adopted to calculate the variables of the ghost cell. When the boundary is moving, a simple temporal extrapolation is used to settle the socalled “freshcell” problem [10]. Test cases show that the method is at least secondorder accurate in space globally and locally.
Moreover, we employ the gaskinetic BGK scheme proposed by Xu [11] to solve the flow field. This numerical scheme is based on the gaskinetic theory and constructed in a more physical way than the NavierStokesfunctionbased scheme. In contrast to the continuum assumption, the physical model of the gaskinetic scheme is more basic and can describe highly nonequilibrium flow, leading to several advantages of this scheme: robust, positivitypreserving and satisfying the entropy condition spontaneously. In smooth flow region the scheme can attain high accuracy while in discontinuous region it can automatically introduce proper viscosity by a physical mechanism so that the discontinuity can be resolved in the resolution of the mesh. Overall, we adopt this scheme because that it is applicable to the viscous flow from incompressible to compressible, especially suitable for the flow involving both strong discontinuity and shear layer, without any complex artificial fixing.
2 Numerical method
In this section, we first sketch the gaskinetic BGK scheme used to govern the flow in our numerical simulation. Then, we describe the grid adopted in our work. Finally, the concept and the construction of our IB method are detailed. It is noteworthy that although the current numerical method is presented in 2D framework, it is not laborious to extend it to 3D situation.
2.1 Gaskinetic BGK scheme
The classical gaskinetic BGK scheme [11] is based on finitevolume framework, in which the whole flow field will be divided into control volumes. As with the traditional finitevolume scheme, the governing equations for the BGK scheme can be written as
(1) 
where is the control volume and is the volume surface. is the vector of conserved variables, defined as
(2) 
where , , and are the mass, momentum and energy densities. is the flux vector. In 2D cartesian grid, Eq. 1 can be discretized as
(3) 
where denotes the average quantity on the cell and is assumed to be located at the cell center to the secondorder spatial accuracy. is the cell area, is the flux at the th interface of the cell, is the length of the interface and is the number of the cell’s interfaces.
The issue of a finitevolume scheme is to calculate the flux at the cell interface. To achieve this point, the scheme with accuracy higher than first order will begin with a reconstruction procedure in which the flow variables around the interface will be constructed. In the classical BGK scheme [11], linear reconstruction is employed and the conserved variable inside the cell is constructed as
(4) 
where and are the slopes of the variable in and directions respectively. Here, in order to adapt the scheme to the flow with shock, the van Leer limiter is used, i.e.
(5) 
where
(6)  
(7)  
(8) 
When implementing this reconstruction in a nonuniform Cartesian grid (detailed later in Subsection 2.2), the procedure is slightly different. Take the grid in Fig. b as an example, for the cell , because the cell is divided into four cells, when calculating by Eq. 7, will be replaced by the average value of the four cells, i.e.
(9) 
By this means, will be a secondorder estimation of the variable slope at center of the left interface of the cell . Besides, for the cell , when calculating the slope at its right interface, due to the absence of the proper adjacent cell, Eq. 6 can not be used. Here, it is calculated by modifying the value of ,
(10) 
where is the size of the cell and can be obtained from slopes at left interfaces of the cell and , i.e.
(11) 
Thus, all of the slopes can be determined in a secondorder manner. In the direction, the slope can be obtained by a same strategy.
After we have constructed the flow around the interface, the flux going across the interface is to be calculated. As a gaskinetic scheme, distinguishing from the traditional NS solver, the BGK scheme takes gas particles as carriers of the conserved quantities. A distribution function is used to describe the particles’ dynamics and the flux is obtained by integrating the particles going across the interface. Let the subscript denote the cell interface of the cell in the positive direction, the relation between the interface flux vector and the distribution function is
(12) 
where and are the particle velocities in and directions respectively. is the internal variable which can be regarded as a vector with degrees of freedom. For 2D flow, the particle motion in direction is included into and is equal to [12]. is the vector of moments
(13) 
is the volume element in the velocity space. For more descriptions about the concepts of these variables please refer to Xu [11].
So the issue is to calculate the gas distribution function at the cell interface. In the BGK scheme, the evolution of the distribution function is governed by the BGK equation [13], whose general solution can be easily obtained. In the classical construction of this scheme, utilizing the Taylor expansion and the ChapmanEnskog expansion, the general solution of the BGK equation is expanded at the interface as [11, 14]
(14)  
Here, different from the previous papers [11, 14], we have rearranged the expression of to make the physical meaning clearer. In Eq. 14, the superscript denotes that the variable is calculated in an upwind manner and will be obtained from variables on different sides of the interface depending on the particle velocity . are Maxwellian distributions and can be got from the conservative variables at both sides of the interface. Other variables, such as , etc, corresponding to the spacial or temporal derivatives, can be obtained from the slopes of the conserved variables at both sides of the interface. More details about these variables can be found in Xu [11] and Xu et al. [14]. Anyhow, surrounding the interface , cells’ information is required in the above calculation and the stencil is shown in Fig. 1.
Since we have obtained the gas distribution function at the interface, the numerical flux across the interface can be computed by Eq. 12, and then the conserved variables can be updated through Eq. 3.
At the end of this subsection, we’d like to give some discussion on the above scheme. In Eq. 14, terms multiplied by represent the initial gas distribution while terms multiplied by represent the nearequilibrium gas distribution. Terms multiplied by are the nonequilibrium parts of the initial and nearequilibrium gas distribution, which are involved with the viscosity and heat conduction. Term multiplied by is the substantial derivative of the equilibrium gas distribution. In the continuum limit, , the gas distribution soon approaches the nearequilibrium state and Eq. 14 can be reduced to
(15) 
which will essentially yield a NS flux [11]. It is good enough to use Eq. 15 to calculate the gas distribution at the interface if the whole flow field is sufficiently smooth. On the other hand, for the flow with discontinuity, the classical BGK scheme introduces an artificial viscosity through
(16) 
where are pressures on two sides of the interface and is the CFL time step. is the physical collision time corresponding to the real physical viscosity while the rest part on the right of the equal sign corresponds to the artificial viscosity. In the flow region which is continuum and smooth, the above artificial viscosity will be negligible and a NS flux (Eq. 15) will be obtained. In the discontinuous region, the collision time will be increased artificially to the scale of and a highly nonequilibrium gas distribution will be obtained through Eq. 14, which means that the thickness of the discontinuous layer will be increased to the scale of the mesh size, leading to the nonoscillating resolving of the discontinuity. This is a welldesigned mechanism based on the basic model of the gaskinetic scheme. The only drawback is that in the case of low cell Reynolds number, the CFL time step is limited, causing the deviation of the gas distribution function from Eq. 15, which will make the scheme break down.
Another thing we want to mention is the Prandtl number fixing. In the BGK model [13], particles with different velocity share a same relaxation rate, resulting in a constant Pr number. Xu [11] directly modified the heat flux to fix the Pr number and a similar idea is adopted here. As stated above, in Eq. 14, terms multiplied by are involved with the heat conduction, hence the heat flux of Eq. 14 can be obtained through a similar way described in Xu [11]
(17)  
where are macroscopic gas velocities obtained from and respectively. Then the energy flux in the vector can be fixed as
(18) 
Here, different from previous works [11], the heat flux is a weighted value calculated from the initial and nearequilibrium gas distributions, which may be more compatible with the general solution Eq. 14.
2.2 Grid strategy
In the IB method, the mesh is not required to fit the solid boundary, thus a nonuniform Cartesian grid is adopted in our work, as shown in Fig. a. To simplify the implementation of the IB method, the grid near the boundary is uniform. A quadtree strategy is used to refine the grid, which means that a cell can be divided into four higherlevel subcells while each of the subcells can be subdivided into four cells with a cell level even higher (see Fig. b). In order to avoid too abrupt grid variation and simplify the data structure, celllevel difference greater than one between neighboring cells is not allowed.
The refinement control can be either artificial or solutionadaptive. In the artificial manner, grid refinement and coarsening are managed by a preset procedure. This is applicable to the problem where the boundary motion is predictable. In this case, grid near the boundary will be refined and the refined region is programmed to move with the boundary. In the solutionadaptive manner, grid refinement and coarsening are controlled by two parameters [15], i.e.
(19) 
(20) 
where is the fluid velocity, is the local cell size, determines the weight of the length scale and is taken to be in the present work. and measure the local strength of the curl and divergence of velocity respectively, which can be used to capture shear layers and shocks [16]. Thus, the refinement criteria is
(21) 
while the coarsening criteria is
(22) 
where , and are userdefined parameters and the value of should be less than to avoid recursive refinement and coarsening.
In a problem where the boundary is stationary and the flow is steady, the above refinement will be done only a few times while in an unsteady problem such procedures will be executed at intervals.
2.3 Immersedboundary method
Forming a mirrored flow field in the solid region to fulfil the boundary condition is applied in plenty of boundaryconforming as well as IB methods [8, 7, 9, 5, 10]. As shown in Fig. 3, in these methods, cells in the flow field are called as “fluid cells”. To complement the calculation stencils of the fluid cells near the boundary, a layer of cells near the boundary and inside the solid body are constructed as “ghost cells” through various of strategies, so that all fluid cells can be updated by a unified scheme.
The same idea is applied in our work. In the present method, the boundary is represented by line segments and a uniform Cartesian grid is used around the boundary. The first step is to identify the fluid cells, whose centers are outside the solid body. The rest of the cells, whose centers lie inside the solid body, are defined as “solid cells”. Although there are various of ways can be used to do this identification including the raycrossing technique [17], here we adopt a way more compatible with our succedent procedure. As shown in Fig. 4, if two cells share a same interface, their cellcenter nodes are connected by dashed line. As a result, a topology of the cellcenter nodes is formed. Then let the solid boundary cut off the connecting lines between different cell centers. After that we can pick an arbitrary fluid cell, all cells whose center nodes are connected, directly or indirectly, with the center of the picked cell are fluid cells. The boundary segments cutting off the cellcenter connecting lines will be recorded and used in the succedent procedure. On the fluid cells, fluid variables will be updated normally by the numerical scheme, while the solid cells will be let alone and not involved in the calculation.
In the next step, the ghost cells will be identified and constructed by a novel strategy. In the previous ghostcell IB methods [8, 9, 5, 10], the ghost cell will be uniquely determined and variables on a ghost cell will have unique values for all interfaces of the cell, which may suffer singularity problem in some cases. As shown in Fig. 5, two boundary segments, and , with different angles, lie in two directions of the ghost cell which can be constructed based on either of the two segments. It is not rational to assume a unique construction for ghost cell in such a case. A reasonable practice is using the information of to construct when calculating the variables at interface , while using the information of to do the construction when dealing with the interface , just as what is done in a conventional bodyfitted grid. Further more, in the conventional bodyfitted ghostcell treatment, when calculating the flux at the interface , the fluid cell which is in the calculation stencil will be regarded as a ghost cell on which the variables will be constructed by the information of and , as if the boundary segment were extended between the cells and . Indeed, the flow passes through the interface will meet the boundary segment , which will exert a reaction on the flow to fulfill the boundary condition. Hence the local boundary segment should determine the construction of the ghost cells when calculating the variables at .
Thus, a “local boundary determination” thought is implemented in this paper to complete the ghostcell construction. When calculating the variable at a interface between a fluid cell and a solid cell, the local boundary segment cutting off the line connecting the two cell centers will be extended to a line. Then this line will be regarded as a virtual solid boundary in the vicinity of the interface. All cells in the calculation stencil whose center nodes lie on the solid side of the line will be identified as ghost cells. Just as what is shown in Fig. 6, when calculating the flux at the interface , the calculation stencil includes cells shown in the picture, where cells, and , are solid cells while the rest cells are fluid cells. The boundary segment which cuts off the connecting line between the cells and will be extended to a virtual boundary line. The solid cells , and the fluid cell whose centers lie on the solid side of the extended will be identified as ghost cells.
After that, the conception of image point [8, 9, 5, 10] is introduced here to construct the variable at the center of the ghost cell. As shown in Fig. 6, for the ghost cell , the point which is symmetric with the cellcenter node about the extended boundary segment is the image point. The line connecting these two symmetric points intersects the extended at the boundaryintercept point , where the boundary conditions for should be satisfied. The boundary condition can be either Dirichlet or Neumann type. For a Dirichlet boundary condition which has the form
(23) 
where is a generic variable at the boundaryintercept point and is the given value for it, the variable at the center node of can be calculated as
(24) 
where is the corresponding variable at the image point . For a Neumann boundary condition, i.e.
(25) 
where is the outer normal vector of , is the gradient at the boundaryintercept point and is the given normal gradient, can be calculated as
(26) 
where is the distance between and the center of . Therefore, once we have obtained the variable at the image point , the variable for the ghost cell can be determined by a linear construction.
Like in the predecessor’s works [5, 10], a bilinear interpolation is adopted to obtain the value of the variable at the image point. The interpolation has the form
(27) 
where are unknown coefficients which can be determined by interpolation conditions. In this step we try to get the real value of the fluid variable at the image point, so the virtual boundary, which is only used for ghost cell determination, is put aside and the information of the real solid boundary is taken into account. Several cases will be encountered as shown in Fig. b. The simplest case is that cells encircling the image point are fluid cells, like what is depicted in Fig. a. In this case, directly substituting the information of the fluid cells into Eq. 27 will yield
(28) 
from which the unknown coefficients can be solved out and the variable at the image point can be determined. The situation will be more complicated when or of the surrounding cells are solid cells. In this case, a line will be drawn between the image point and the center of the solid cell, which will intersect the boundary segment at a point where the boundary condition can be used to close the interpolation. Just as what is shown in Fig. b, the line segments between the image point and the center nodes of the solid cells intersect the boundary segments and at and respectively. If a Dirichlet boundary condition is applied at and , a function similar to Eq. 28 can be got where the variables corresponding to the solid cells and will be replaced by the corresponding values at and . If a Neumann boundary condition is applied, i.e.
(29)  
where and are the outer normal vectors of and , taking the derivative of Eq. 27 and substituting it into the boundary conditions Eq. 29, and then combining the interpolation conditions at the fluid cells will yield
(30) 
where and are coordinates of and respectively. Afterwards the unknown coefficients can be solved out and the value of the variable at the image point can be obtained by the interpolation Eq. 27.
Thus, just like what is shown in Fig. 6, all ghost cells in the calculation stencil can be determined and eventually all fluid cells can be updated by a unified scheme.
When implementing the above method in a simulation where the boundary is moving, the socalled “freshcell” problem [10] will be encountered. This happens when a solid cell turning into a fluid cell along with the boundary motion. Just as displayed in Fig. 8, the cell is a solid cell at the th, th steps and turns into a fluid cell at the th step. How to determine the variable on this new fluid cell (i.e. fresh cell) is the issue. Here a simple extrapolation in time is adopted to resolve this problem. A generic variable for the fresh cell at th step will be calculated as
(31) 
where is the time and the superscript denotes the number of the time step, and are the ghostcell values of the cell at th, th steps respectively. It is notable that in the present method, a solid cell may corresponds to more than one ghost cell, whose variables will be constructed from different boundary segments. So and in Eq. 31 are obtained in a weighting manner, i.e. if the cell at the th step has ghostcell constructions corresponding to the boundary segments respectively, then can be calculated as
(32) 
where is the normal distance that the center of outside the boundary segment at the th step, is the ghostcell value corresponding to , is the Heaviside function defined as
(33) 
The weight for certain ghostcell value is calculated in consideration of the fact that the center node of the fresh cell will be biased to a more distant (in the sense of normal distance) boundary segment, based on which the ghostcell value should have a larger weight. As shown in Fig. a, the center node of the cell is biased to which is more distant than (i.e. ) at the th step, thereby a larger weight will be given to the ghostcell value constructed from . In Fig. b, the center node of is totally biased to and lies on the solid side of at the th step, which means that the center node never goes over to the fluid side of and is always a point inside the solid body for this boundary segment. In this case a negative will be got and the weight for the ghostcell value based on is .
Another notable thing is that, just like various of previous IB methods [5, 10], the present IB method cannot guarantee the conservative condition near the immersed boundary, which may lead to a persistent increase or decrease in the conserved quantity. Therefore it is necessary to give boundary conditions which can make the problem wellposed. If the flow is wellposed, the conservative problem may become just a matter of accuracy.
3 Numerical results and discussions
In this section, several test cases are conducted to validate the present method. First, the spacial accuracy of the method is investigated in the simulation of TaylorCouette flow. Then, the test case of supersonic flow over a stationary circular cylinder is carried out to testify the ability of the present method in handling highly compressible flow. Finally, the flow over an oscillating circular cylinder and the flow over a moving airfoil are conducted to demonstrate the capability of the method in managing the moving boundary.
In this part, the force which fluid exerts on the solid body is calculated by summing the momentum fluxes going from the fluid cells into the solid cells. Just as shown in Fig. 10, the direction component of the force can be calculated as
(34) 
where is the time step, and are the momentum and mass fluxes going across the th interface into the solid cell, is the velocity of the solid body and is the area of the interface.
3.1 TaylorCouette flow
We have simulated the TaylorCouette flow to assess the spatial accuracy of the current method. In this test case, viscous fluid is filled in the gap of two concentric rotating cylinders, as shown in Fig. 11. In incompressible limit, an analytical solution can be solved out for this flow [18], and the variables at radius can be calculated as
(35) 
(36) 
(37) 
where is the circumferential velocity of the fluid, is the pressure, is the ratio of specific heat, is a variable involved in the gaskinetic scheme and related to the gas temperature which can be calculated as , where is the specific gas constant. The integration constants can be determined by the radiuses and the boundary conditions of the inner and outer cylinders, which are omitted here for simplicity.
Our simulation employs a square computational domain with uniform Cartesian grid, as shown in Fig. 11. The radiuses of the two cylinders are set to be . The inner cylinder is holding stationary while the outer cylinder is rotating at a velocity equal to (approximating the incompressible flow). The Reynolds number defined by the width of the gap and the velocity is . The Prandtl number is set to be . To make the problem wellposed, the boundary conditions at the two cylinders are set to be
(38) 
where all of the variables are nondimensionalized by the density and the speed of sound at the outer cylinder. grid sizes, , are adopted to simulate this flow.
First, the classical BGK scheme with a full expansion of the distribution function Eq .14 is employed to solve the flow field. The and error norms of the velocity , pressure and (related to the temperature) comparing with the analytical solution for different spatial resolutions are shown in Fig. 12, in which the calculation of the error norms involves all of the fluid cells in the flow field. It is demonstrated that, for the velocity, the measure accuracy is around and measure accuracy is , both are beyond the expected secondorder accuracy, although in our method all interpolations and constructions are no more than secondorder spacial accuracy. This may result from the socalled “superconvergence” of the classical BGK scheme [19]. It’s worth pointing out that although the convergence rate is faster, for the coarse grid the error of BGK scheme seems larger than that of the conventional NS solver. For the pressure, the spacial accuracy is in sense while only in sense, and the curve for error norm is somewhat zigzag. The pressure errors on the fluid cells near the inner cylinder are shown in Fig. a. The figure demonstrates that very high pressure error will occur on the fluid cell part of which lies on the solid side of the immersed boundary, which implies that the data from these cells may be not reliable. Another notable thing is that this problem is much lighter near the outer cylinder where the boundary is concave. The maximum pressure error of the fluid cell near the outer cylinder is one order of magnitude smaller than that near the inner cylinder. Thus, we think these singular fluid cells are due to the singularity of the “ghost flow field”, i.e. one piece of flow field lies on the solid side of the immersed boundary (referred to as the “ghost flow field”) may corresponding to several pieces of real flow field lie on the fluid side thus variables on the ghost flow field may have singular values, which is not suffered by a concave boundary. This is just similar to the singularity of the ghost cell discussed in Subsection 2.3. However, these singular fluid cells seem not to spoil the order of the local accuracy as we exclude the fluid cells adjoining one or more solid cells from the calculation of the error norms and show the result in Fig. 13. It can be seen that, the accuracy of pressure achieves in sense and in sense. For , as shown in Fig. 12 or Fig. 13, the spacial accuracies in and measure are both straightly firstorder. Furthermore, the error contours (Fig. b) show that the max error occurs near the middle of the gap between the two cylinders, which implies that it is the BGK scheme but not the IB method which gives rise to the low accuracy of .
Then, we reduce the full expansion of the distribution function (Eq .14) to Eq. 15 and construct a simplified BGK solver to do some further investigation. It is a fact that a linear interpolation will give us a secondorder estimation of the variable itself while a firstorder estimation about the slope of the variable. Therefore, for Eq. 15, if the linear reconstruction is employed just as what is done in the classical BGK scheme, the slopeassociated terms, i.e. the terms multiplied by , will be firstorder accurate in space. In general, in the continuum limit, and this firstorder error will be covered up. However, when we use this simplified BGK scheme with linear reconstruction to do the above simulation of TaylorCouette flow, we get a result similar to what is shown in Fig. 13 and firstorder accuracy of is acquired. After that, we add compensatory terms to Eq. 15 to fix the firstorder errors of the terms multiplied by and use this modified BGK scheme to do the same simulation. The and error norms for different spacial resolutions are shown in Fig. 15. It can be seen that, for the accuracy of velocity or pressure, there is no significant difference from the previous result (Fig. 13), but the spacial accuracy of is increased to in sense and in sense. We have made an analysis for this and put forward the idea that it is the specificity of this test case which causes the low accuracy when using the classical BGK scheme with linear reconstruction. In the TaylorCouette flow, the heat convection is ignorable and the temperature of the fluid, which is related to , is dominated by the viscosity and heat conduction, which are only determined by the slopeassociated terms (i.e. the nonequilibrium part of the gas distribution function) in the fullexpansion distribution function Eq .14 or the simplified distribution function Eq. 15. Thus, if one wants to get secondorder accuracy about the temperature (or ) in such a heatconvectionignorable test case, the slopes of the variables should be reconstructed in a secondorder manner.
Finally, for the IB method, we can conclude that the local spacial accuracy (in sense, excluding the fluid cells adjacent to solid cells) is secondorder and the global spacial accuracy (in sense) is also at least secondorder.
3.2 Supersonic flow around a stationary circular cylinder
The supersonic turbulent flow past a stationary circular cylinder has been considered as a test case to show the performance of the present method on simulating highly compressible flow. The freestream Mach number is . The Reynolds number based on the freestream conditions and the cylinder diameter is . The Prandtl number is set to be . A rectangular computational domain with a size of is adopted where the cylinder is placed at . In the far field, a coarse grid with the size is used. Approaching the cylinder, the grid is gradually refined and the finest grid is used around the immersed boundary. Besides, the solutionadaptive local grid refinement is employed and the grid will be refined at shear layers and shocks at every time step, as shown in Fig. 16. The flow field is directly solved by the classical gaskinetic BGK scheme without any turbulencesimulation technique, which makes the simulation very challenging. In view of the computational cost, three sizes of the nearwall grid, , and are considered, although even the finest size seems not refined enough to resolve the turbulent flow.
The instaneous flow pattern and the corresponding vorticity field, pressure coefficient contours and Mach number contours obtained from the finest grid are shown in Fig. b, Fig. 18 and Fig. 19 respectively. Despite of the insufficient spatial resolution, some details of the flow can be observed. As what is displayed in the figures, a bow shock is formed ahead of the cylinder and there is a subsonic region behind the shock. The subsonic stream expands and accelerates along the surface of the cylinder, soon developing to the supersonic flow, separating at the rear of the the cylinder. Then a subsonic turbulent recirculation region forms behind the cylinder, enveloped by the separated supersonic flow. A strong shear layer can be observed between the subsonic recirculation region and the supersonic flow (Fig. b). Finally the supersonic streams from the upper and lower surfaces of the cylinder converge, forming two oblique shocks in the wake. The Mach number contours (Fig. 19) look very compatible with the published data in [4] and [6].
The mean vorticity distributions nondimensionalized by the freestream velocity and the diameter of the cylinder for different nearwall grid sizes are shown in Fig. 20. It can be seen that the vorticity curves are suffered from some oscillations. These oscillations are mainly due to the inadequate spatial resolution near the solid boundary and the way we obtain the vorticity, i.e. directly calculating the velocity gradient at the boundaryintercept point from the velocity of the corresponding image point, in which the velocity of the image point may be constructed oscillatory along the solid boundary if the grid is not refined enough. Despite of the oscillations, it can be observed that the vorticity increases significantly as the grid is refined, which implies that the solution may be still far from convergence at the finest grid. Besides, the locations of separation for different nearwall spatial resolutions determined from Fig. 20 are given in Tab. 1, where the total drag coefficients are also listed. The present results are compared with the experimental result of Ignatova and Karimullin [20] and the numerical results of De Palma et al. [4] and de Tullio et al. [6]. It can be found that for the drag coefficient , the present results are consistent with the experimental result and also agree reasonably well with other numerical results. For the location of separation , the present results appear nonconvergent and deviate from most of the experiment and numerical results by around except the numerical result obtained by De Palma et al. [4] from their coarse grid, which agrees well with the present .
Work  
Present ()  1.422  
Present ()  1.436  
Present ()  1.446  
Ignatova and Karimullin [20]  1.43^{*}  
De Palma et al. (IB method, coarse) [4]  1.38  
De Palma et al. (body fitted, fine) [4]  1.40  
De Tullio et al. [6]  1.41  

The experimental pressure drag coefficient in [20] is taken as the total drag coefficient approximately.
The mean pressure coefficient distributions for different nearwall spatial resolutions are plotted and compared with other authors’ results in Fig. b. Fig. a shows that as the grid is refined, the oscillation of the pressure coefficient distribution is suppressed. This oscillation about the pressure can be also seen in the result of De Palma et al. [4] obtained by their IB method, which seems like a common problem for the class of IB methods. Besides, it can be observed that as the spatial resolution increases, the pressure at the rear of the cylinder will decrease, leading to the increase in total drag (Tab. 1). In Fig. b, the present result and the numerical result of De Palma et al. [4] obtained from their bodyfitted mesh almost overlap, showing great consistence, and they are also in good agreement with the numerical result of de Tullio et al. [6]. All three numerical results agree with the experimental data reasonably well. The most significant difference between the experimental and numerical results is that at up to the separation point the experimental curve lies below the numerical curves. This difference can be attributed to the spatial effect in the experiment [20], which is ignored in the twodimensional numerical simulation.
In general, for this test case, without any turbulencesimulation technique and in the insufficient spatial resolution, the present method fails to give a smooth vorticity distribution and a very accurate separation point. However, there are still a lot of flow characteristics can be resolved correctly. With proper turbulence model or LES technique, we are confident that the present method can yield a satisfactory result for this supersonic turbulent problem.
3.3 Incompressible flow around an oscillating circular cylinder
The test case of a circular cylinder oscillating incompressibly in stationary fluid has been conducted to demonstrate the ability of the present method in manipulating the moving boundary. In this test case, a circular cylinder oscillates sinusoidally in the horizontal direction with the equation of motion which can be expressed as
(39) 
where is the horizontal displacement of the cylinder, is the amplitude and is the frequency. Two key parameters dominating the pattern of this flow are the Reynolds number Re and the KeuleganCarpenter number KC, which are defined as
(40)  
where is the diameter of the cylinder, is the maximum horizontal velocity of the cylinder, is the fluid density and is the fluid viscosity. In our simulation, the condition is considered and the maximum velocity is set to be equal to a Mach number to approximate the incompressible flow. According to Tatsuno and Bearman [21] and our previous work [22], for this parameter set the flow is twodimensional, stable and symmetric. When the cylinder moves through the equilibrium position, the maximum velocity is attained and two counterrotating vortices with the same magnitude of strength are formed behind the cylinder. Then the cylinder will reach the maximum displacement and move back, hitting at the pair of vortices which will be split and finally reverse during the backward motion of the cylinder.
As shown in Fig. 22, a rectangular computational domain with the size is adopted in this test case. The equilibrium position of the cylinder’s sinusoidal motion is set at . A nonuniform Cartesian grid, as described in Subsection 2.2, is adopted with a large mesh size equal to in the far field and a refined size (refined enough [22]) near the cylinder. The refined region around the immersed boundary will move along with the motion of the cylinder, as what is displayed in Fig. 23. The classical gaskinetic BGK scheme is employed in the simulation.
The instantaneous flow pattern at phase angle is shown in Fig. 24. At this moment the cylinder has just reached the maximum negative displacement and is moving back to its equilibrium position, hitting at the two vortices formed behind it. It can be seen that the stream line fits the boundary completely and no flow penetration can be observed.
The instantaneous pressure and vorticity isolines at phase positions and are shown in Fig. 25 and Fig. 26 respectively. In the left figures the present result is compared quantitatively with the result from Yuan et al. [22] (our previous work), where isolines from two results share the identical variable values. As is demonstrated in the left figures, the two results coincide with each other perfectly near the cylinder while minor mismatches of the isolines exist in the area a little further. These mismatches are attributable to the insufficient computing time in Yuan et al. [22], which will make the flow nonfully developed in the far field. In the right figures of Fig. 25 and Fig. 26, the numerical result of Dütsch et al. [23] is exhibited (variable values of the isolines may be different from those in the left figures) and we can get some qualitative comparison between these three results. Generally speaking, the three results agree reasonably well.




The velocity profiles near the cylinder at four vertical cross sections for three phase positions , and are shown in Fig. 27. Here, and are defined as
(41) 
where are the coordinates relative to the equilibrium position of the sinusoidal oscillation, are the velocity components in the horizontal and vertical direction. The present result is compared with the numerical and experimental results of Dütsch et al. [23]. Generally good agreement is obtained among the three sets of results.



According to Morison et al.’s semiempirical equation [24], which is widely used in estimating the inline force on the body in oscillatory flow, the timedependent inline force on the cylinder oscillating in a stationary fluid can be expressed as
(42) 
where is the displacement of the cylinder, and are the drag coefficient and the added mass coefficient. and are two empirical hydrodynamic coefficients and can be obtained by various methods such as Fourier analysis or leastsquares fitting after we have got the timedependent inline force from experiment or numerical simulation. In our numerical simulation, the inline force is obtained from Eq. 34 and the values of and fitted to this force are shown and compared with other authors’ numerical results in Tab. 2. The present result agrees well with other three sets of results. The inline force calculated from the Morison equation Eq. 42 with the present and is shown in Fig. 28 and compared with the forces obtained directly from the numerical simulations of us and Dütsch et al. [23]. It can be seen that the force calculated from the Morison equation is grossly consistent with the two numerical results, which agree with each other excellently.
3.4 Compressible flow around a moving airfoil
To further verify the ability of the present method in handling compressible moving boundary problem, the test case of an airfoil passing through a compressible viscous fluid is conducted. In this simulation, a NACA0012 airfoil at an angle of attack is moving at a Mach number in a stationary fluid. The Reynolds number based on the chord length and the freestream condition is . The airfoil is set to be insulated and the Prandtl number is . A rectangle computational domain with a large size is adopted, which will move along with the motion of the airfoil so that the airfoil will be always placed at the center of the computational domain. In the far field a coarse grid with a size of is used while near the airfoil the grid is gradually refined to , as shown in Fig. 29. Also, the refined region will move with the airfoil. The classical gaskinetic BGK scheme is employed in the simulation. Furthermore, as a control, the flow around a stationary NACA0012 airfoil with the same conditions, computational domain and grid is investigated.
Fig. a shows the stream lines around the airfoil in the reference frame which is fixed in the space and the reference frame which is fixed on the airfoil. In Fig. a, it can be seen that the two sets of flow patterns, which are based on the movingairfoil simulation and the stationaryairfoil simulation respectively, are entirely coincident and fit the solid boundary perfectly. It is demonstrated that the flow separates at some place from the upper surface of the airfoil and a large recirculation region forms behind the airfoil. The location of the separation , which is determined from the distribution of the skin friction coefficient (Fig. 34), the drag coefficient and the lift coefficient , for which the drag and lift forces can be obtained by Eq. 34, are presented in Tab. 3. It shows that the two sets of present results agree completely and they are also consistent with the results from other authors.
Work  
Present (moving)  0.2734  0.4347  0.363 
Present (stationary)  0.2734  0.4347  0.363 
Müller et al. [26]  0.2632  0.4199  0.371 
Cambier [27]  0.2656  0.4342  0.36 
Kordulla [28]  0.2845  0.4261  0.362 
Paillere and Deconinck [29]  0.2728  0.4530  0.383 
The Mach number contours and the pressure coefficient contours from the movingairfoil simulation and the stationaryairfoil simulation are shown and compared with other numerical results in Fig. 31 and Fig. 32. Excellent agreement is acquired between the two sets of present results, and they seem also coincident with the numerical results of Müller et al. [26] and Hafez and Guo (Overflow) [30].
Due to the imposition of the adiabatic boundary condition, the temperature will vary along the surface of the airfoil. The distributions of based on the moving airfoil and the stationary airfoil are plotted in Fig. 33, where passible agreement is obtained between the two curves. It can be observed that at the leading edge the temperature is around or even higher than the free stream total temperature , although it is lower than at most surface of the airfoil.