obstacle interactions using an improved ghost-cell immersed boundary method

Computers and Fluids 182 (2019) 128–143

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

Benchmark solutions

Numerical simulations of shock/obstacle interactions using an improved ghost-cell immersed boundary method Yang Zhang a, Xinglong Fang b, Jianfeng Zou a,∗, Xing Shi a, Zhenhai Ma a, Yao Zheng a a b

Center for Engineering and Scientific Computation, and School of Aeronautics and Astronautics, Zhejiang University, Hangzhou 310027, China AECC Hunan Aviation Powerplant Research Institute, Zhuzhou 412002, China

a r t i c l e

i n f o

Article history: Received 29 August 2018 Revised 4 December 2018 Accepted 14 February 2019 Available online 15 February 2019 Keywords: shock/obstacle interaction boundary layer immersed boundary method

a b s t r a c t An improved ghost-cell immersed boundary method coupled with a high-order finite difference scheme is currently proposed to numerically simulate the reflected shock/boundary layer interaction in a wavy-wall tube. The main improvement of current algorithm involves the particular treatment of the solid boundary reconstruction procedure by using both ghost points in solid domain and forcing points in fluid domain together. Owing to this kind of improvement operation, the amount of points used to recover the solid boundary increases greatly and will inevitably promote the accuracy of solid boundary reconstruction, which will play a significant role in the current high-order accurate numerical simulation. A typical validation case of 2-dimensional shock/circular cylinder interaction is conducted and demonstrate that the current improved numerical method can be used to capture the flow details such as shock bifurcation and small vortex very well and improve the numerical precision greatly. The improved method is then applied to numerically investigate the effects of boundary layer on the propagation of the shock in a wavy-wall tube. It is observed that, at low Mach number, the existence of boundary layer has little effect on the shock propagation. As the Mach number increases to 1.9, some certain complicated flow features will occur due to the enhanced boundary layer effects. In the early stage, the flow structures are nearly the same as those without the boundary layers. In the transition stage, a high-speed shear band behind the leading shock is observed and finally evolves into a worm-like structure, whose instability is illustrated by qualitative assessment on its flow features. In the developing stage, a remarkable shock bifurcation phenomenon is present and the worm-like structures will start shedding regularly, which will induce a series of scattering-like waves and asymmetrical flow fields. These characteristics demonstrate that the bifurcated shock waves located in upper and lower walls of the wavy-wall shock tube act as a primary trigger to gradually destabilize the shocked flow. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction The immersed boundary method (IBM) was firstly introduced by Peskin [1] to model the blood flow in the heart. The main idea of this proposed method is to introduce a virtual forcing term in order to represent the interaction between the heart walls and the fluid on uniform Cartesian grids. More recently, Goldstein et al. [2] and Saiki and Biringen [3] proposed a feedback forcing model to accurately represent the effect of solid body. However, this kind of ‘forcing term’ approaches will inevitably bring spurious numerical oscillations and the time step is severely limited if an explicit time-stepping scheme is adopted. In addition, a penalty method was introduced by Arquis and Caltagirone [4] and a few improved versions were published [5,6]. These approaches add ∗

Corresponding author. E-mail address: [email protected] (J. Zou).

https://doi.org/10.1016/j.compfluid.2019.02.014 0045-7930/© 2019 Elsevier Ltd. All rights reserved.

additional forces to the continuous N-S equation. In contrast, Mohd-Yusof [7] firstly proposed a direct forcing method for the simulations of problems with rigid boundaries. This kind of forcing is directly considered in the context of the discretized equation to drive the numerical solution towards the required boundary conditions and no additional terms are added in the governing equations, which will successfully solve the stiffness problem [8]. Fadlun et al. [9] combined this approach with a finite-difference method on a staggered grid and showed that this method was more efficient than feedback forcing model. The basic idea of ghost-cell immersed boundary method is similar to that proposed by Ghias et al. [10]. It characterizes the boundary geometry information indirectly through ghost points within the solid geometry to satisfy the corresponding boundary conditions. And there is also no need to add additional items to the flow governing equations. Therefore, the ghost-cell IBM can be easily combined with the existing in-house CFD solver, and presents

Y. Zhang, X. Fang and J. Zou et al. / Computers and Fluids 182 (2019) 128–143

wide application in dealing with flows with complex geometry. Mittal et al. [11] designed a sharp interface immersed boundary method based on a ghost-cell scheme to solve three-dimensional incompressible flows of complex geometry. In this newly proposed method, the image point related to each ghost cell is achieved by symmetrically mirroring it across the boundary, and a bilinear interpolation from the surrounding fluid points available is applied to identify the flow state at the image point. This proposed IBM method is claimed to be second-order accurate in space compared with the previous cut-cell interface methods [12,13]. To further improve the accuracy of the boundary treatment, in the high-order immersed-boundary method currently proposed by Zhu et al. [14], the novel interpolation procedure was performed directly with the first layer of fluid points and a third-order polynomial was used to ensure at least fourth-order accuracy of the flow reconstruction at the forcing point. Actually, the most critical issue concerned for ghost-cell IBM is the reconstruction procedure, which is used to accurately recover the flow information of the so-called ghost points on geometric boundary. Usually, the flow states of the image points along the normal direction of the boundary are reconstructed by a certain interpolation algorithm from the surrounding fluid points [10,11]. The Neumann boundary conditions are frequently used to avoid negative weighting coefficients in reconstruction process [15]. Tseng and Ferziger [16] systematically investigated the advantages and disadvantages of first-order linear reconstruction, quadratic polynomial reconstruction and inverse distance weighted reconstruction, respectively, in dealing with either Dirichlet or Neumann boundary conditions. Although the higher-order reconstruction of ghost points is able to preserve higher resolution of flow simulation, it will always increase the numerical instability [17]. In addition to reconstruction scheme, the way to identify ghost-cell nodes is also of great importance on the features of convergence and accuracy for the ghost-cell IBM. The relatively large stencil of fluid cells will always lead to much more ghost-cells embedded deeply inside the immersed body [10]. Chaudhuri et al. [18] applied 2-layer ghost-cells to match the Weighted Essentially Non-Oscillatory (WENO) scheme, in order to obtain enough interpolation points for solving the shock/obstacle interactions. It’s apparent that the more the layers of ghost-cells close to the solid boundary are defined, the more the geometric boundary is identified with high resolution. Besides, mesh refinement near the boundary is also required to meet the demand of large stencils for high-order scheme near the boundary, which will make it possible to achieve accurate flow solution near the boundary. However, it should be noted that both mesh refinement operation and increasing ghost-cell layers will inevitably lead to a more expensive computation, especially for the shock/obstacle interaction problems involving complex shock/boundary layer and shock/shock interactions. Qu et al. [17] proposed a kind of constrained moving leastsquares (MLS) reconstruction method to recover the flow states of the ghost points and successfully eliminates the numerical instability encountered in high-order simulation for the interaction of supersonic viscous fluids (Ma = 1.2–3.0) with arbitrarily shaped moving solids. Bailoor et al. [19] developed a fluid-structure interaction (FSI) solver for simulating the impact of a blast wave on thin elastic structures. A high-order polynomial interpolation combined with a weighted-least square error minimization proposed by Seo and Mittal [20] is used to impose boundary conditions. A feedback penalty immersed boundary method combined with a fifthorder weighted essentially non-oscillatory (WENO) scheme was adopted to solve the fluid-structure interaction with compressible multiphase flows involving large structure deformations and discontinuities of shock waves are clearly captured in this work [21]. Luo et al. [22] tested the advantages of ghost-cell boundary

129

method in calculating compressible flows, with varying Reynolds numbers (<50 0 0) and Mach numbers (<0.5), respectively. A ghostpoint-forcing/direct-forcing hybrid method proposed by Boukharfane et al. [23] has a capability to resolve a large variety of complicated flows, such as subsonic laminar viscous flows and high-speed shock-obstacle interactions. Owing to the great benefit from the use of the Cartesian grid system in the context of ghost-cell IBM, it is very easy to apply high-order finite-difference schemes and implement massively parallel programming. Hence the ghost-cell IBM is becoming the most promising numerical method for solving supersonic or hypersonic compressible flow problems with complicated configurations in geometry, which was used in the present numerical investigation on the shock/obstacle interactions. In this paper, the author proposed an improved ghost-cell IBM based on the work of Chaudhuri et al. [18], which is coupled with a high-order finite difference scheme to achieve temporally accurate numerical solution of the supersonic compressible flow problems. In the present study, two distinct numerical cases of shock/obstacle interactions were calculated to demonstrate the advantages of the currently improved method in capturing the shock and vortex structures. The first case is about the 2-dimensional shock/cylinder interaction. The second case is the study of the interaction between the moving reflected shock and the boundary layer in a wavy-wall tube. The most primary feature of the second validated problem is the shock bifurcation phenomenon. It was firstly discovered by Mark [24] when studying the propagation characteristics of shock waves in a end-closed pipe. Laterly, Davies and Wilson [25], Matsuo et al. [26] and Fokeev [27] studied the geometric characteristics of the shock bifurcation structure in details. Due to the high speed of the reflected shock wave, the shock/boundary layer interaction is a remarkable dynamic process and presents strong unsteadiness. Weber et al. [28] pointed out that there exists a well marked instability problem in the shock/boundary layer interaction. Gamezo et al. [29,30] analyzed the impacts of the moving shock bifurcation on flame acceleration. Penyazkov and Skilandz [31] experimentally studied the effects of wall roughness on the characteristics of the shock bifurcation. They found that small-scale turbulent structures caused by small roughness is able to stabilize the shock structure. The present research was a further extension of the previous work performed by Lodato et al. [32]. In Londato’s study, the simulation of a plane shock interacting with a sinusoidal wavy-wall (wall amplitude and wavelength are 1 mm and 20 mm, respectively) was conducted. They analyzed the dynamic evolution of the triple point of λ-shock, the transverse wave system and the critical points in details. And the unstable propagation characteristics of the reflected shock wave were also discussed. Different from the previous studies by Lodato et al. [32], a non-slip adiabatic condition instead of a periodic boundary was imposed on the side walls of the tube in our present simulations. Thus the propagation of the reflected shock front and the transverse waves are strictly confined inside the closed tube, which could be a more realistic scenario in engineering. The shock bifurcation is then induced near the side walls owing to the reflected shock/boundary layer interaction, which will certainly increase the instability and complexity of the shocked flows. The primary purpose of the current research is to reveal the specific features of the evolution of the shock bifurcation in the tube with the presence of wavy ends and non-slip adiabatic walls, which is of great significance to understand the intrinsic physical mechanism related to the shocked flow. The organization structure of this paper is given as follows: an improved ghost-cell immersed boundary method is proposed in Section 2. A 2-dimensional shock/circular cylinder interaction (Ma = 2.81) is taken as a validation case in Section 3 to test the accuracy of the proposed method. The interactions between the reflected shock

130

Y. Zhang, X. Fang and J. Zou et al. / Computers and Fluids 182 (2019) 128–143

point (IP). The information on the ghost points is then easily determined by a second-order linear approximation, which will effectively avoid occurs of negative values for weighting coefficient.

φG =2φB −φI

Fig. 1. Basic immersed boundary method.

wave and the boundary layer in a 3D wavy-wall tube are simulated and analyzed in details in Section 4. Finally, a summary is provided in Section 5. 2. Mathematical calculation method 2.1. Improved ghost-cell immersed boundary method According to Chaudhuri et al. [18], as a numerical preparation, it is of essential role to identify the ghost nodes within the solid domain and store the necessary parameters related to these nodes, such as element properties, geometric coordinates and flow states, prior to the start of a new simulation. As shown in Fig. 1, the entire computational domain used currently consists of a solid domain (solid ) and a fluid domain (fluid ). To distinguish between the solid domain and the fluid domain, a set of marked variables (ϕ i, j ) is defined on all Cartesian grid points (xi , yj ).

fluid = solid =





ϕi, j = 1, (xi , y j ) ∈ fluid , i = 1, ..., Nx, j = 1, ..., Ny



 ϕi, j = 0, (xi , y j ) ∈ solid , i = 1, ..., Nx, j = 1, ..., Ny



GP =

(xi , y j ) ∈ solid if ∃(xk , yl ) ∈ fluid

φF =(φB + φI )/2

φI =

4  m=1

ω m φm 

ω m = ηm

(3)

The information of ghost points (φ G ) should be reconstructed in a certain manner through the fluid nodes to satisfy the requirement of the imposed physical boundary conditions. The imaging approach is applied in the current simulation to define the image

(5)

It’s usual that the location of the image point (IP) does not coincide exactly with the mesh point, so the information of IP needs to be reconstructed from the surrounding fluid points. The most frequently used method to fulfill the reconstruction process is the linear interpolation. Luo et al. [22] takes the four fluid nodes of the cell where IP is located inside and use bilinear interpolation to obtain the flow information of IP. Majumdar et al. [15] reported that there is not remarkable difference in precision between bilinear interpolation and quadratic interpolation. Chaudhuri et al. reconstructed the information of IP by an inverse distance weighted interpolation method proposed by Franke [33], which ensures the smoothness of the reconstruction solution near the extreme point.

(2)



for k = i − 3, ..., i + 3, l = j − 3, ..., j + 3

In practice, the number of grid points around the geometric boundary will certainly decide the amount of image points and ghost points achieved, and the accuracy of boundary identification as well. Local mesh refinement is then required to retrieve highresolved boundary geometry and accurate boundary flows. However, it should be noted that both mesh refinement operation and increasing ghost-cell layers will lead to a more expensive computation. In order to solve this problem, a compromise is reached between the high resolution and the expensive computation. Currently, an additional layer of forcing points (FP) close to the boundary in fluid domain is marked, and same linear interpolation of the mirror relationship is used to reconstruct these forcing points. Fig. 2 makes a comparison about the identification of the boundary curve AB with or without this additional layer of forcing points. Due to the existence of this layer of forcing points, the number of curve segments for boundary identification increases to 15 (nearly two times of that of the case with only ghost points), and the boundary point (BP) distribution tends to be more uniform than the previous method with only ghost points.

(1)

Herein Nx and Ny are the number of grid points in x and y directions, and fluid , solid represent the point sets of fluid domain and embedded solid domain, respectively. To successfully resolve the shock discontinuity and the complicated turbulent details in the supersonic flows, a fifth-order WENO and sixth-order central difference schemes are adopted to model the convection term and the viscous term in Navier–Stokes equations, respectively. To meet the requirement of the number of stencil points for the accurate high-order discretization near the boundary, three layers of ghost-cell points are identified around the solid boundary. The ghost points (GP) can be marked as follows:

(4)

4 

(6)

 ηk

(7)

k=1

ηm = 1/dk2

(8)

where φm (m = 1, 2, 3, 4 ) denote physical information on the four neighbor points (NP) in the cell containing IP, and dk is the distance from the kth NP to IP, as shown in Fig. 1. It should be noted that a particular case will arise from the use of additional forcing point, i.e. the degenerate case in which the forcing point FP coincides with one of its four neighbor points NP. In this particular case, the number of neighbor points NP collected to perform interpolation for the forcing point is reduced to three. In order to work around this kind of degenerate situation, the four-point inverse distance weighted interpolation method is extended to a multi-point interpolation within a prescribed distance. As shown in Fig. 3, the fluid points trapped in a circle centered at the image point IP with a radius of a prescribed value R are identified as the neighbor points to recover the flow states of the forcing or ghost points.

φI =

np  m=1

ωm φm (np is the number of neighbor points )

(9)

Y. Zhang, X. Fang and J. Zou et al. / Computers and Fluids 182 (2019) 128–143

131

Fig. 2. Subdivision of solid boundary.

series) and 12 cores. The memory of a single node is up to 64GB to ensure a high computational efficiency. 3. Shock/cylinder interaction

Fig. 3. Improved immersed boundary method.

2.2. High-order finite difference solver The current version of immersed boundary method does not require a specially customized fluid solver. Due to the use of Cartesian grids, current flow simulations are performed by a high-order finite difference method coupled with the improved ghost-cell IBM introduced above. The Steger-Warming flux vector splitting method [34] combined with a fifth-order accurate weighted essentially non-oscillatory (WENO) scheme is used to deal with the convective terms. A sixth-order accurate central difference scheme for the diffusion terms and a three-step Runge-Kutta approach for time advancing are adopted. The numerical accuracy of this DNS solver used in our work has been widely demonstrated in the investigation of the supersonic flows [35,36]. Current numerical calculations are achieved in the National Supercomputing Center of Tian He in Guangzhou City in China. In this high-performance computing platform, each computer node is equipped with two CPUs (Xeon E5

As the supersonic flow problems are concerned, the shock reflection has been widely studied in the previous research [37–39]. When an incident shock wave encounters a cylindrical obstacle during propagation, a reflected shock occurs around the surface of the cylinder to form a converging 2-shock structure. As the interaction point keeps moving downstream over the surface of the cylinder, this 2-shock structure will soon transform into a typical Mach reflection, i.e. a λ-shock bifurcated structure with a triple point. The shock/cylinder interaction in 2D space is currently used to validate the current ghost-cell immersed boundary method. The plane incident shock Mach number is 2.81, the length and width of the whole calculation domain are Lx × Ly = 240mm × 60mm respectively, and the diameter of the cylinder D is 10 mm. The supersonic inflow and outflow are imposed to the left and right boundary respectively, and the upper and lower walls are both nonslip adiabatic walls. The entire calculation domain adopts uniform Cartesian grids, and totally five cases with different grid size are conducted to confirm the grid independence, as shown in Table 1. Fig. 4 shows the contours at different time instants to illustrate the evolution of the shock/cylinder interaction. The primary flow features, i.e. the incident shock (I.S.), the reflected shock (RS), the Mach stem (M.S.1 and M.S.2), the triple points (T.P.1 and T.P.2), and the slip lines (S.L.1 and S.L.2), are clearly observed. A weak vortex structure (V.) can also be observed near the slip lines. Especially in the wake behind the cylinder, the upper and lower Mach stems collide with each other, and a secondary λ-shock structure arises Table 1 Grid size with 5 different cases. Case

Nx × Ny

x = y

Mesh1 Mesh2 Mesh3 Mesh4 Mesh5

500 × 125 667 × 167 10 0 0 × 250 20 0 0 × 500 40 0 0 × 10 0 0

0.048D 0.036D 0.024D 0.012D 0.006D

132

Y. Zhang, X. Fang and J. Zou et al. / Computers and Fluids 182 (2019) 128–143

Fig. 5. Pressure coefficients of the cylindrical surface in different grid sizes.

Fig. 4. Schlieren images for the shock reflection evolution (upper: numerical results, bottom: experimental data [40]).

in the wake (see Fig. 4d). Compared with the experimental visualization of Bryson and Gross [40], our currently simulated results successfully capture the characteristics of the flow field very well. Fig. 5 shows the pressure coefficients of four locations on the cylindrical surface (ρ 0 and u0 is the supersonic inlet density and velocity, respectively). The number ‘1’ in this figure denotes the stagnation point of the front surface, and ‘2’, ‘3’, ‘4’ are three points with 60°, 120°, 180° measured in a clockwise direction from the front stagnation point, respectively. It is known that as the mesh size is refined, the numerical solution gradually converges, and the mesh independence meets the resolution requirement in case of Mesh4 and Mesh5. In the spatial resolution of Mesh1, there is a small numerical oscillation at the stagnation point 1. In the flow field behind the cylinder, due to the shock/wake and shock/shock interactions, the flow structures become extremely complicated and failed to be highly resolved in Mesh1 and Mesh2. Due to the high resolution, the flow solution in 0.03 ≤ t/t0 ≤ 0.16 for Mesh5 can be used as a benchmark for the unsteady shock/cylinder interaction problem. The L2 norm of the local error is expressed as

ε=

1 Nt

0.16  t/t0 =0.03





benchmark 2 t

ψt − ψ

1/2 (10)

Fig. 6. L2 norm of the error for pressure coefficient versus the grid size.

where Nt represents the number of time steps in t/t0 ∈ [0.03, 0.16], and Nt = 25 is adopted here.ψ t is the instantaneous physical quantity of the flow field. In Fig. 6 we present a log-log plot of the errors of pressure coefficients for point 2 and 3 on the cylindrical surface as a function of the grid size. It is observed that these error norms show a convergence accuracy close to second-order, that is to say the current improved immersed boundary solver is locally second-order accurate in terms of solid boundary treatment. As mentioned previously, the currently proposed ghost-cell immersed boundary method is mainly improved by considering the ghost point (GP) in the solid domain and the forcing point (FP) in the fluid domain together to recover a high-resolved boundary geometry. Table 2 lists six different cases, where GP method and GP+FP method mean the boundary condition is exerted by

Table 2 Calculation conditions in different IBMs. Case

Nx × Ny

Methods

C1 C2 C3 C4 C5 C6

667 × 167 10 0 0 × 250 20 0 0 × 500 667 × 167 10 0 0 × 250 20 0 0 × 500

GP+FP method GP+FP method GP+FP method GP method GP method GP method

Y. Zhang, X. Fang and J. Zou et al. / Computers and Fluids 182 (2019) 128–143

Fig. 7. Pressure coefficients with different IBMs.

ghost points only and ghost/forcing points together, respectively. The case of C3 is selected as a benchmark according to the test of grid independence. Fig. 7 shows a comparison of pressure coefficients achieved by the currently improved ghost-cell method with those by the original ghost-cell method. It can be seen that for the improved ghostcell method namely GP+FP method, the relative numerical error is distinctly smaller than the original ghost-cell method (GP method). For both cases of C1 and C4 (coarse grids), the pressure coefficient at the stagnation point of the cylindrical front shows little difference from the reference solution. But as the point-3 and point4 are concerned, there is a remarkable deviation from the reference solution for both two methods. For the cases of C2 and C5, the improved ghost-cell method managed to achieve a higher resolution to the boundary flow around the cylinder compared with that of the original ghost-cell method. As the number of grids increases further (see C3 and C6), the resolved flow field does not improve significantly, which demonstrates the mesh independence is approached. Fig. 8 presents the trajectories of the triple points (TP1 and TP2) for the cases of C2, C3, and C5, respectively. The calculated results of the reference case (C3) are in good agreements with the experimental data from Bryson and Gross [40] and Kaca [41], which illustrates that the current high-order finite difference

Fig. 8. Comparison of triple point trajectory.

133

134

Y. Zhang, X. Fang and J. Zou et al. / Computers and Fluids 182 (2019) 128–143

Fig. 9. Streamline diagrams (left: internal treatment of solid; right: no treatment).

based improved ghost-cell boundary method can be used to accurately simulate the shock/obstacle interaction. Additionally, Fig. 9 is the streamline distribution around the cylinder for the case of C3 at t = 22.9 μs. The points of the solid domain inside the cylinder are imposed with a zero speed (left in Fig. 9) or treated in the way that is exactly the same as those points in the fluid domain (right in Fig. 9). It can be observed that the calculation results of the fluid flow are same with each other between these two different treating methods for the internal solid nodes, which illustrates that it’s not necessary to put any particular efforts to deal with the points inside the solid domain, and agrees well with the claims proposed by Iaccarino and Verzicco [8]. 4. Reflected shock/boundary layer interaction 4.1. Computational domain and initial condition According to the experimental data [32,42], the entire computational domain is a three-dimensional channel with a square section of width × height of 2 cm × 2 cm and with a length of 16 cm. Fig. 10 depicts a geometrical diagram of its x-y section. The left boundary in the x-direction is a no-slip adiabatic solid wall of a shape of sinusoidal curve with the amplitude and wavelength of 1 mm and 20 mm, respectively. Uniform Cartesian grids are covered in the whole domain and the wavy-wall of the left boundary is fully embedded, which enable the use of the improved ghostcell immersed boundary method to achieve high-resolution solution with such an easily generated but high-quality grid. In the work of [32], the grid resolution is of 222.2–266.7 μm, which is not fine enough to resolve the flow structures clearly, and a non-

physical contact discontinuity due to numerical error was observed in their simulation. Hence, a more refined grid size of 100 μm is applied in present simulation to ensure the spatial resolution of DNS. A plane shock wave propagating from right to left is initially located at x = 12.88 cm. The right boundary is a supersonic inflow of air, and the region between the right boundary and the shock wave is patched with the uniform post-shock fluid. Given the incident shock Mach number (M1 ), the post-shock parameters can be determined according to the Rankine-Hugoniot relationships. Regarding the boundary conditions, the front and back boundaries in the span-wise z direction are periodic. In order to study the effects of the wall boundary layer on the shock structure, totally four cases corresponding to different incident shock Mach numbers and boundary conditions in the normal y direction are performed, as shown in Table 3. Case 1 is selected as a benchmark and its specific physical parameters are shown in Table 4. Note that superscript ∗ indicates dimensionless for all variables in the following data analysis.

Table 3 Four calculation conditions. Case

Incident shock Mach

y-direction boundary condition

1 2 3 4

1.5 1.5 1.9 1.9

Period No-slip wall Period No-slip wall

Table 4 Flow initial parameters.

Fig. 10. Computation domain schematic.

Quantity

Symbol

Value

Domain size Number of elements Wavy-wall amplitude Wavy-wall wavelength Incident shock Mach Incident shock speed Initial speed Initial density Initial pressure Viscosity (at 291.15 K)

Lx × Ly × Lz Nx × Ny × Nz Aww

16.0 × 2.0 × 2.0 (cm) 1600 × 200 × 200 1.0 mm 2.0 cm 1.5 −514.0 m/s 0.0 m/s −237.96 m/s 1.208 kg/m3 2.25 kg/m3 101,325.0Pa 249,090.6Pa 1.827 × 10−5 kg/(m·s)

λww M1 ui u1 u2

ρ1 ρ2

p1 p2

μ0

Y. Zhang, X. Fang and J. Zou et al. / Computers and Fluids 182 (2019) 128–143

135

Fig. 11. Contours of normalized density after reflection (IS, incident shock; MS, Mach stem; RS, reflected shock; SL, slip line).

4.2. Numerical verification in the benchmark Firstly, the basic flow field is analyzed for the benchmark of C1. As shown in Fig. 11, the plane shock wave is propagating towards the left boundary and reflected by the wavy-wall, resulting in a reflected shock wave with a smooth curve. There is no Mach stem and triple point observed initially (such as the time instant of 6 μs). This is due to the small amplitude of the curved wall, which acts only as a shallow reflector. However, at the time instant of 30 μs, a triple point appears on the front of the reflected shock clearly and a typical λ-shock of Mach reflection is also observed, which consists of a leading shock wave (IS), a Mach stem (MS) and a reflected shock wave (RS). A series of reflected shock waves propagating backwards in the shocked flow will form a crossed wave system (in [32] called weak transverse waves), and a slip line (SL) behind the triple point extends backwards in the shock gas. As the leading shock wave gradually advances towards the right end, the triple point undergoes an unsteady transverse motion. It will cause a series of transverse waves and slip lines, which will interact with each other in the shock gas and form a cellular pattern similar to the detonation flow (such as 102 μs). Fig. 12 gives a qualitative comparison of the simulated density contours with the experimental schlierens. The results are presented at the three instants of 120 μs, 200 μs and 280 μs, respectively, after the incident shock wave is reflected by the wavy-wall. It can be seen that the propagation of transverse shock waves in the shock gas and the motion of the triple points agree well with the experimental observations. At the initial moment, the propagation speed ui of incident shock wave is 514 m/s, and the time required to reach the left wall is about 250 μs. Subsequently it will take 280 μs for the shock wave to approach the location shown in Fig. 12c. At this moment, the distance between the leading shock and the wavy-wall is about 9.1 cm, and the average propagation speed of the leading shock wave (ur ) is 325 m/s. According to the state of the upstream gas with (u, ρ , p) = (−237.96 m/s, 2.25 kg/m3, 249,090.6 Pa), it can be calculated that the Mach num-

Fig. 12. numerical results and experimental data [41].

ber of the reflected shock wave at this moment is 1.43. In addition, by measuring the displacement change of the triple point in the y direction on the leading shock front, the average transverse speed of the triple point is up to 310 m/s. These quantitative results are all in good agreements with the experimental data [42]. 4.3. Weak boundary layer effect (M1 = 1.5) 4.3.1. Basic flow features in Case 1 and Case 2 Fig. 13 shows the density contours for Case 1 and Case 2 with periodic and non-slip adiabatic boundary conditions in y direction, respectively. The leading shock wave (IS) will interact with the boundary layer (see Fig. 13a), which is resulted from the presence of the solid wall. It can be observed that, for the case of C2, the leading shock is curved with a very small curvature (marked by A and B in Fig. 13) near the top and bottom walls where the

136

Y. Zhang, X. Fang and J. Zou et al. / Computers and Fluids 182 (2019) 128–143

Fig. 15. Flow velocity distribution normalized by wall shear velocity before reflected shock.

Fig. 13. Contours of normalized density after reflection (upper: Case 1; bottom: Case 2).

density changes slightly. But the crossed wave system formed by the impacting between transverse waves and walls in the shock gas is basically the same for Case 1 and Case 2. This indicates that the wall boundary layer has little effect on the entire flow field when the incident shock intensity M1 is 1.5. In the following, the flow structure near the wall where the shock wave interacts with the boundary layer in Case 3 is analyzed in details. 4.3.2. Specific flow features near the wall Fig. 14 shows the distribution of the dimensionless velocity u∗ at the time instant of 270 μs in the section of z = 1 cm. There is a high-speed region where the stream-wise velocity is positive near the wall in the shock gas, which is called the “high-speed shear band”. The transverse waves propagate in the shock gas and collide with the high-speed shear band near the wall, which will form a crossed wave system. On the right side of the reflected shock wave, the flow state is same to the supersonic inflow to propagate in the tube. Fig. 15 shows its velocity distribution normalized by the wall shear velocity in three different locations. According to the boundary layer properties measured at x = 10 cm, the spatial resolution in wall units is y+ w = 1.2. The stream-wise velocity increases linearly along the wall-normal direction, which indicates a laminar flow before the reflected shock/boundary layer in-

Fig. 16. Flow velocity distributions of shear bands at different locations.

teraction occurs. This is because the inflow has not fully developed into turbulence in such a short time. The velocity distributions in the high-speed shear band in four different locations are described in Fig. 16. These symmetrical curves indicate the almost identical shear bands for both upper and lower walls. The stream-wise velocity gradually increases along the x direction from x = 4.8 cm to 8 cm, and the locations for the maximum in the upper and lower shear bands are at y = 0.0134 cm and 1.9866 cm, respectively. The stream-wise velocity decreases rapidly from a positive maximum peak to a negative minimum value in a short distance, thereby forming a region with large velocity gradient in the shear band. In addition, the normal velocity v∗ in the shear bands is close to 0. As shown in Fig. 17, the curve denotes the distribution of velocity gradient is at the position of x = 7.27 cm. The velocity gradient outside the shear band is almost zero, and approaches the peak

Fig. 14. Contour map of stream-wise velocity at 270 μs.

Y. Zhang, X. Fang and J. Zou et al. / Computers and Fluids 182 (2019) 128–143

137

Fig. 19. Evolution of maximal velocity gradient with time in the shear band.

Fig. 17. Curve of stream-wise velocity gradient in the lower shear band.

Fig. 20. Curve of shear layer thickness at 270 μs. Fig. 18. Variation of maximal velocity gradient in the shear band. Table 5 Parameter changes of power function. t(μs)

c

d

190 230 270

0.033 0.0013 0.0 0 03

3.636 5.051 5.468

It has been observed in Fig. 14 that the upper and lower shear bands are the extremely thin areas close to the boundary walls, and the flow velocity changes very sharply in these regions. In order to quantify the deformation in the shape of the shear band, a quantity about the thickness δ of the shear band is defined as:

δt (x )= ∂ u∗ /∂ y∗

value = 11 at y = 0.0517 cm. Due to the existence of the transverse shock, the curve presents a slight bump at y = 0.56 cm. Fig. 18 shows the variation of the peak of velocity gradient in the high-speed shear band along the x direction. It can be seen that the closer to the leading shock wave, the larger the velocity gradient peak. And it exhibits an incremental law of power function with the x-direction displacement, as follows:

∂ u∗ ∂ y∗



= c xd ( c > 0, d > 0 )

(11)

umax − umin

(∂ u/∂ y )max

(12)

umax and umin are the maximum and minimum stream-wise velocity in the normal y-direction, respectively. (∂ u/∂ y)max is the peak value of the velocity gradient at the position of x. Fig. 20 depicts the variation of the thickness of the lower shear band at t = 270 μs. It can be seen that the thickness of the entire shear band fluctuates within the range of 0.014–0.02 cm, which implies the inherent instability in the high-speed shear bands. Hence, the boundary layer effect is weak and the instability effect is not remarkable yet at M1 = 1.5.

max

The constant coefficients c and d will change with time and the larger the value of d, the faster the curve grows. Fig. 19 shows the growing of the maximal velocity gradient in the shear bands at t = 190 μs, 230 μs, 270 μs. As the time is advancing, the slope of the curve is increasing. It means that the velocity variation in the shear bands becomes more and more complicated as the leading shock wave gradually pushes forward. As shown in Table 5, the values of the constant coefficients c and d are obtained by fitting these curves. At 270 μs, the constant coefficient d reaches its largest value.

4.4. Strong boundary layer effect (M1 = 1.9) 4.4.1. Basic flow features in Case 3 and Case 4 Fig. 21 shows the temporal evolution of the density contours for Case 3 and Case 4. For the no-slip wall boundary condition (right), the entire flow fields become more complicated. At the early stages (such as 27.89 μs), the shock front (IS and MS) and the subsequent transverse wave (RS) are approximately the same as the left case (Case 3), except that the high-speed shear bands near the upper and lower walls (marked by L). The thickness of the current shear

138

Y. Zhang, X. Fang and J. Zou et al. / Computers and Fluids 182 (2019) 128–143

Fig. 21. Comparison of density contours (left: Case 3; right: Case 4).

band is distinctly larger compared with that in M1 = 1.5. At the time of 47.26 μs, the slip line (SL) collides with the shear band near the wall and multiple slip lines occur due to complex wave reflections. These slip lines interact with the reflected transverse wave (RS), resulting in a bending break in the middle of the RS (see marked circle at 47.26 μs). The motion of the high-speed shear band is followed by the leading shock wave moving to the right. On the one hand, the bending degree of the shock front near the wall is gradually increasing (such as IS at 81.14 μs). And on the other hand, the high-speed shear band slowly evolves into a recirculation zone with a small number of vortices shedding from the wall (marked by V at 81.14 μs). At that time, although the reflected transverse wave RS on the shock front collides with the separated vortices, the crossed waves and the slip lines in the shock gas are basically the same as those in the periodic boundary condition. As the leading shock continues to propagate towards the right boundary, a shock bifurcation appears in the shock front near the wall (such as B in 119.86 μs), and a large number of vortex structures in the recirculation zone are sequentially detached from the walls and shed to the central region. The shock gas is fully filled with many scattering-like structures (such as 119.86 μs and 148.89 μs), which will interact with the crossed transverse waves and cause an extremely non-uniform flow structure. As the Case 4 of the non-slip wall boundary condition is concerned, the interaction process between the leading shock wave and the boundary layer is mainly divided into three stages: the early stage (before the time of 47.26 μs), the transitional stage (47.26 μs–81.14 μs) and the development stage (after 119.86 μs), respectively. At this early stage, the shock front is less curved near the wall. In the vicinity of the wall, there is a high-speed shear band. In the transitional stage (47.26 μs–81.14 μs), the high-speed shear band gradually evolves into a recirculation zone with a small amount of vortex shedding, and the transverse waves in the shock gas are clearly observed, and are basically the same as those of the case with periodic boundary condition. In the development stage

(after 119.86 μs), the leading shock wave is bifurcated near the wall, and the vortex continuously detaches from the wall. Scattered wave system appears in the shock gas, and the flow field is extremely complicated. 4.4.2. Early stage The overall flow field in the early stage is similar to that at M1 = 1.5, but there are still many differences. Fig. 22 depicts the dimensionless velocity distributions along the normal direction at three locations (x = 0.95, 1.13, 1.33 cm) in the high-speed shear bands at 47.26 μs [u∗ = (−0.28, 1.3), v∗ = (−0.22, 0.22)]. The maximum stream-wise velocity u∗ is greater than the inflow velocity of 1.0, and the average magnitude of fluctuation is greater than that at M1 = 1.5 [u∗ = (−0.12, 0.33), v∗ = (−0.038, 0.038)]. A curve in Fig. 23 shows a change of the thickness of the lower shear band that increases first and then decreases. The maximum and the minimum are 0.007 cm to 0.03 cm, respectively. These results indicate a more complicated shear band when the Mach number of the incident shock increases. Figs. 24 and 25 show the trajectory of triple point and the shape of leading shock front, respectively. The trajectory of triple point or the front shape is closely coincident with that of the case with periodic boundary condition, except that for the flow near the wall which is affected by the boundary layer effect. As shown in Fig. 21, in the early stage, the crossed transverse waves are basically the same in both cases. These indicate that at the early stage, the interaction of the wall boundary layer with the shock wave does not significantly affect the transverse propagation of the triple point on the shock front. Based on the above analysis, when the initial incident shock Mach number is increased, the interaction between the leading shock wave and the boundary layer is enhanced a lot. However, in the early stage, the flow structures (including the shape of the shear band, the velocity field distribution, and the fluctuation of the thickness of the high shear band) are all close

Y. Zhang, X. Fang and J. Zou et al. / Computers and Fluids 182 (2019) 128–143

139

Fig. 22. Flow velocity distributions of the shear band at 47.26 μs.

Fig. 23. Curve of shear band thickness at 47.26 μs.

Fig. 25. Leading shock front shape.

Fig. 24. Trajectory of trip-point. Fig. 26. Contour maps of stream-wise velocity at 47.26 μs, 61.78 μs, 76.3 μs, 95.66 μs.

to those of the low Mach number case. The transverse motion of the triple point on the front shock caused by the left wavy-wall boundary is not significantly affected by the boundary layer.

4.4.3. Transitional stage In the transitional stage, the high-speed shear band gradually evolves into a recirculation zone, and a small amount of vortex fall off from the walls. As shown in Fig. 26, the streamwise velocity

140

Y. Zhang, X. Fang and J. Zou et al. / Computers and Fluids 182 (2019) 128–143

Fig. 27. Shape change of shock front in the transition stage.

distributions at four moments clearly illustrate that the shape of the high-speed shear band gradually changes from a smooth curve to a worm-like (w) geometry. This worm-like structure is filled with the positive high-speed fluid and the negative reverse flow is observed in the region just around this worm structure to meet the requirement of dynamic balance (such as circle mark at 76.3 μs). The maximum reverse velocity is up to u∗ = −0.495, which is close to half of the speed of the inflow. As the leading shock wave propagates forwards, the bending degree of the shock front near the wall gradually increases (shown in Fig. 27). The shape of the shock front is transformed from a smooth curve to several connected line segments, which is identified as the shock bifurcation (marked by B).

As shown in Fig. 28, the flow velocity u∗ distribution is zoomed in at 76.3 μs, from which the worm-like high-speed shear band can be clearly observed. The shock front is a smooth curve and only a pair of vortices is detached from the wall in the shock gas. These structures will be specifically analyzed below. As shown in Fig. 29, the velocity distributions in the high-speed shear band at different streamwise positions are almost the same as those in the early stage. The streamwise velocity rapidly decreases after a large peak value in a short distance, then another local peak point is reached, and finally the flow velocity in the central region approaches zero. At different x positions, the maximum value or the minimum value and the corresponding positions of these extreme values obtained in the normal y-direction are all different with each other. At x = 2.03 cm, the maximum and minimum values are obtained at y = 0.027 cm and y = 0.124 cm, respectively. At x = 2.17 cm, the extreme points are obtained at y = 0.075 cm and y = 0.198 cm, respectively. The difference of the normal locations of the extreme points for different streamwise positions will indicate the distortion of the internal structure of the high-speed shear band. Figs. 30 and 31 are the local friction coefficient and pressure distribution on the lower wall at 66.6 μs and 76.3 μs. It can be seen that in the upstream region away from the worm-like structure, the wall friction coefficient is small and close to zero. When it is approaching to the high-speed ‘worm’ structure (at 66.6 μs corresponding position x = 1.5 cm) that the wall friction coefficient suddenly increases and then stats to oscillate. At 66.6 μs, the maximum amplitude of the fluctuation is about 0.47C¯ f and there is an abrupt decrease for the friction coefficient as the leading shock is passing through. The wall pressure also drops abruptly at the upstream part of the worm-like structure and irregular fluctuations occur inside the shear band. These unsteady fluctuations of the wall pressure and friction coefficient in the high-speed shear band demonstrate the inherent instability of the curved worm-like structure.

Fig. 28. Enlarged view of stream-wise velocity at 76.3 μs.

Fig. 29. Flow velocity distributions of the shear band at 76.3 μs.

Y. Zhang, X. Fang and J. Zou et al. / Computers and Fluids 182 (2019) 128–143

Fig. 30. Friction coefficient distribution.

141

Fig. 33. Leading shock front in the development stage.

Fig. 31. Static pressure distribution.

4.4.4. Development stage The primary features in the development stage are described as: (1) Shock front near the wall is bifurcated; (2) Continuous vortex shedding phenomenon occurs in the tail of the recirculation region; (3) Many scattering-like wave structures are present in the shock gas and interact with the crossed transverse waves. Fig. 32 displays the flow contours of two physical quantities at the time instant of 119.86 μs. As can be seen from the pressure distribution, the entire leading shock front can be roughly divided into three parts: the leading shock bifurcated foot (IS1), the curved segment shock (IS2) and the straight Mach stem (MS). The curved segment shock wave is the transitional section of the shock front and connects with the upper and lower triples. The upper triple point is created by the shock bifurcation in the wall boundary layer, and the lower triple point is created by the collision of the incident shock wave with the left end wall due to the presence of the left wavy-wall. As the leading shock wave is advancing,

the upper triple point gradually moves toward the center line, and the height of the bifurcated foot increases slowly. The lower triple point moves up and down along the shock front, which will result in a crossed transverse wave system in the shock gas. Fig. 33 shows the evolution of the shock wave’s shape in the periodic and no-slip wall boundary conditions together. The circle mark in Fig. 33 represents the triple point of the transverse wave, and the triangle mark is the shock bifurcation point near the wall. It can be found that the trajectory of the triple point in the transverse direction is nearly the same for both cases. This implies that the boundary layer flow does not significantly affect the transverse motion of the triple point on the shock front. However, the propagation speed of the case with no-slip wall boundary will gradually differ from that of the case with periodic boundary as the time is advancing, which is mainly impacted by the flow field behind the shock wave. The shock bifurcation is ready to accelerate the propagation of the leading shock wave in the presence of the boundary layer. It can also be seen from the velocity gradient distribution in Fig. 32 that behind the λ-shock bifurcation, there exists a continuous vortex shedding from the tail of the recirculation zone. Similar vortex shedding process was also previously observed in the studies by Bulovich et al. [43] and Daru and Tenaud [44]. Many scattering-like wave structures ( marked in the Fig. 32) appear in the shock gas, which are resulted from the vortex shedding process, and interact with the reflected transverse shocks (RS). The whole flow field is no longer symmetrical due to the presence of the shock bifurcation and vortex shedding processes. 5. Conclusion This paper proposes an improved ghost-cell immersed boundary method, which is combined with a high-order finite difference algorithm, to numerically simulate the supersonic flows with

Fig. 32. Contours of flow variables at 119.86 μs.

142

Y. Zhang, X. Fang and J. Zou et al. / Computers and Fluids 182 (2019) 128–143

complex geometry. The way to treat the solid boundary in the present improved approach differs from that of the original ghostcell method. Currently, the ghost points inside the solid domain and the forcing points inside the fluid domain are applied together to accurately reconstruct the effective boundary around the solid domain. As the verification case of a 2D shock reflection problem is concerned, currently proposed approach was successfully used to achieve a higher resolution of the shocked flow, compared with the general ghost-cell IBM proposed by Chaudhuri et al. [18]. The shock wave/boundary layer interactions in a reflected wavy-wall are numerically studied, where the moving shock wave and the transverse wave system are strictly confined in the tube. It is revealed that when the incident shock Mach number is 1.5, the shock/boundary layer interaction is weak, and the wall boundary layer has a little effect on the main flows. When the incident shock Mach number increases to 1.9, the shock/boundary layer interaction is enhanced significantly. Based on the primary features and the dynamics of the flow field, the whole process is divided into three stages: early stage, transition stage and development stage. In the early stage, the shock front is less curved near the wall, and the flows in the high-speed shear band are close to those in the low Mach number case. In the transition stage, the highspeed shear band gradually evolves into a worm-like recirculation zone. The measured normal locations of the local peak values of the velocity gradient for different streamwise positions differ from each other. The observed fluctuations of the shear band thickness, wall pressure and friction coefficient clearly indicate the inherent instability of such worm-like structure. However, the transverse waves in the shock gas are still visible, and are nearly the same as those in the case of periodic boundary condition. In the development stage, the trajectory of the triple point is basically consistent with that in the periodic boundary condition. This indicates that the role of the boundary layer does not significantly affect the transverse motion of the triple point of the shock front. The shock bifurcation arises near the wall, and there is a continuous vortex shedding phenomenon at the tail of the worm-like structure. The shedding vortex firmly interacts with the transverse shock waves, which will induce scattering-like wave structures and a greatly inhomogeneous and asymmetrical flow field. Acknowledgments The National Supercomputing Center of Tianhe in Guangzhou provides efficient computing resources for simulations in this paper. Financial support from the National Natural Science Foundation of China (Grant no. 11372276, 11432013, 11372278) is sincerely acknowledged. References [1] Peskin CS. Flow patterns around heart valves: a numerical method. J Comput Phys 1972;10(2):252–71. [2] Goldstein D, Handler R, Sirovich L. Modeling a no-slip flow boundary with an external force field. J Comput Phys 1993;105(2):354–66. [3] Saiki EM, Biringen S. Numerical simulation of a cylinder in uniform flow: application of a virtual boundary method. J Comput Phys 1996;123(2):450–65. [4] Arquis E, Caltagirone JP. Sur les conditions hydrodynamiques au voisinage d’une interface milieu fluide-milieu poreux: applicationa la convection naturelle. CR Acad. Sci. Paris II 1984;299:1–4. [5] Kolomenskiy D, Schneider K. A Fourier spectral method for the Navier–Stokes equations with volume penalization for moving solid obstacles. J Comput Phys 2009;228(16):5687–709. [6] Gazzola M, Chatelain P, Van Rees WM, Koumoutsakos P. Simulations of single and multiple swimmers with non-divergence free deforming geometries. J Comput Phys 2011;230(19):7093–114. [7] Mohd-Yosuf J. Combined immersed boundary/ B-spline methods for simulations of flow in complex geometries. Annual Research Briefs, Center for Turbulence Research; 1997. p. 317–28.

[8] Iaccarino G, Verzicco R. Immersed boundary technique for turbulent flow simulations. Appl Mech Rev 2003;56(3):331–47. [9] Fadlun EA, Verzicco R, Orlandi P, Mohd-Yusof J. Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations. J Comput Phys 20 0 0;161(1):35–60. [10] Ghias R, Mittal R, Dong H. A sharp interface immersed boundary method for compressible viscous flows. J Comput Phys 2007;225(1):528–53. [11] Mittal R, Dong H, Bozkurttas M, Najjar FM, Vargas A, von Loebbecke A. A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries. J Comput Phys 2008;227(10):4825–52. [12] Udaykumar HS, Shyy W, Rao MM. Elafint: a mixed Eulerian-Lagrangian method for fluid flows with complex and moving boundaries. Int J Numer Methods Fluids 1996;22(8):691–712. [13] Udaykumar HS, Mittal R, Rampunggoon P, Khanna A. A sharp interface Cartesian grid method for simulating flows with complex moving boundaries. J Comput Phys 2001;174(1):345–80. [14] Zhu C, Luo H, Li G. High-order immersed-boundary method for incompressible flows. AIAA J 2016:2734–41. [15] Majumdar S, Iaccarino G, Durbin P. RANS solvers with adaptive structured boundary non-conforming grids. Annual Research Briefs, Center for Turbulence Research, Stanford University; 2001. p. 353–466. [16] Tseng YH, Ferziger JH. A ghost-cell immersed boundary method for flow in complex geometry. J Comput Phys 2003;192(2):593–623. [17] Qu Y, Shi R, Batra RC. An immersed boundary formulation for simulating high-speed compressible viscous flows with moving solids. J Comput Phys 2018;354:672–91. [18] Chaudhuri A, Hadjadj A, Chinnayya A. On the use of immersed boundary methods for shock/obstacle interactions. J Comput Phys 2011;230(5):1731–48. [19] Bailoor S, Annangi A, Seo JH, Bhardwaj R. Fluid–structure interaction solver for compressible flows with applications to blast loading on thin elastic structures. Appl Math Modell 2017;52:470–92. [20] Seo JH, Mittal R. A high-order immersed boundary method for acoustic wave scattering and low-Mach number flow-induced sound in complex geometries. J Comput Phys 2011;230(4):10 0 0–19. [21] Wang L, Currao GM, Han F, Neely AJ, Young J, Tian FB. An immersed boundary method for fluid–structure interaction with compressible multiphase flows. J Comput Phys 2017;346:131–51. [22] Luo K, Zhuang Z, Fan J, Haugen NEL. A ghost-cell immersed boundary method for simulations of heat transfer in compressible flows under different boundary conditions. Int J Heat Mass Transfer 2016;92:708–17. [23] Boukharfane R, Ribeiro FHE, Bouali Z, Mura A. A combined ghost-point-forcing/direct-forcing immersed boundary method (IBM) for compressible flow simulations. Comput Fluids 2018;162:91–112. [24] Mark H. The interaction of a reflected shock wave with the boundary layer in a shock tube. National Advisory Committee for Aeronautics; 1958. [25] Davies L, Wilson JL. Influence of Reflected Shock and Boundary-Layer Interaction on Shock-Tube Flows. Phys Fluids 1969;12(5):1–37. [26] Matsuo K, Kawagoe S, Kage K. The interaction of a reflected shock wave with the boundary layer in a shock tube. Bull JSME 1974;17(110):1039–46. [27] Fokeev VP. Peculiarities of reflected shock wave bifurcation on a needle. Tech Phys Lett 2015;41(8):793–6. [28] Weber YS, Oran ES, Boris JP, Anderson JD Jr. The numerical simulation of shock bifurcation near the end wall of a shock tube. Phys Fluids 1995;7(10):2475–88. [29] Gamezo VN, Khokhlov AM, Oran ES. The influence of shock bifurcations on shock-flame in-teractions and DDT. Combust Flame 2001;126(4):1810–26. [30] Gamezo VN, Oran ES, Khokhlov AM. Three-dimensional reactive shock bifurcations. Proc Combust Inst 2005;30(2):1841–7. [31] Penyazkov O, Skilandz A. Bifurcation parameters of a reflected shock wave in cylindrical channels of different roughnesses. Shock Waves 2018;28(2):299–309. [32] Lodato G, Vervisch L, Clavin P. Direct numerical simulation of shock wavy-wall interaction: analysis of cellular shock structures and flow patterns. J Fluid Mech 2016;789:221–58. [33] Franke R. Scattered data interpolation: tests of some methods. Math Comput 1982;38(157):181–200. [34] Steger JL, Warming RF. Flux vector splitting of the inviscid gasdynamic equations with application to finite-difference methods. J Comput Phys 1981;40(2):263–93. [35] Li X, Fu D, Ma Y. Direct numerical simulation of hypersonic boundary layer transition over a blunt cone with a small angle of attack. Phys Fluids 2010;22(2):025105. [36] Zhang Y, Zou J, Xie J, Li X, Ma Z, Zheng Y. Asymmetric characteristics of the shock bifurcation in the reflected shock/boundary layer interaction. Int J Numer Methods Heat Fluid Flow 2018;28(10):2357–77. [37] Giepman RHM, Schrijer FFJ, Van Oudheusden BW. Flow control of an oblique shock wave reflection with micro-ramp vortex generators: effects of location and size. Phys Fluids 2014;26(6):066101. [38] Fang J, Yao Y, Zheltovodov AA, Lu L. Investigation of three-dimensional shock wave/turbulent-boundary-layer interaction initiated by a single fin. AIAA J 2016;55(2):509–23. [39] Igra D, Takayama K, Igra O. Shock wave reflection over roughened wedges. In: 30th international symposium on shock waves, 1. Springer; 2017. p. 621–5. [40] Bryson AE, Gross RWF. Diffraction of strong shocks by cones, cylinders, and spheres. J Fluid Mech 1961;10(1):1–16.

Y. Zhang, X. Fang and J. Zou et al. / Computers and Fluids 182 (2019) 128–143 [41] Kaca J. An interferometric investigation of the diffraction of a planar shock wave over a semicircular cylinder; 1988. NASA STI/Recon Technical Report N, 89. [42] Briscoe MG, Kovitz AA. Experimental and theoretical study of the stability of plane shock waves reflected normally from perturbed flat walls. J Fluid Mech 1968;31(3):529–46.

143

[43] Bulovich SV, Vikola˘ınen VÉ, Zverintsev SV, Petrov RL. Numerical simulation of the interaction between reflected shock wave and near-wall boundary layer. Tech Phys Lett 2007;33(2):173–5. [44] Daru V, Tenaud C. Numerical simulation of the viscous shock tube problem by using a high resolution monotonicity-preserving scheme. Comput Fluids 2009;38(3):664–76.