obstacle interactions

Journal of Computational Physics 230 (2011) 1731–1748

Contents lists available at ScienceDirect

Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp

On the use of immersed boundary methods for shock/obstacle interactions Arnab Chaudhuri, Abdellah Hadjadj ⇑, Ashwin Chinnayya CORIA UMR 6614 CNRS – INSA, University of Rouen, 76801 Saint-Etienne du Rouvray, France

a r t i c l e

i n f o

Article history: Received 4 March 2010 Received in revised form 14 July 2010 Accepted 9 November 2010 Available online 2 December 2010 Keywords: Immersed boundary method Shock–obstacle interaction WENO scheme

a b s t r a c t This paper describes the implementation of immersed boundary method using the directforcing concept to investigate complex shock–obstacle interactions. An interpolation algorithm is developed for more stable boundary conditions with easier implementation procedure. The values of the fluid variables at the embedded ghost-cells are obtained using a local quadratic scheme which involves the neighboring fluid nodes. Detailed discussions of the method are presented on the interpolation of flow variables, direct-forcing of ghost cells, resolution of immersed-boundary points and internal treatment. The method is then applied to a high-order WENO scheme to simulate the complex fluid–solid interactions. The developed solver is first validated against the theoretical solutions of supersonic flow past triangular prism and circular cylinder. Simulated results for test cases with moving shocks are further compared with the previous experimental results of literature in terms of triple-point trajectory and vortex evolution. Excellent agreement is obtained showing the accuracy and the capability of the proposed method for solving complex strongshock/obstacle interactions for both stationary and moving shock waves. Ó 2010 Elsevier Inc. All rights reserved.

1. Introduction Studies of shock-wave interaction involving complex geometries is an important field of research not only in fundamental sciences but also for engineering applications. For instance, the interaction of shocks with solid obstacles can occur in several applications including high-speed flows with particles, hypervelocity impact and penetration. These phenomena have received significant attention in the past and yet remain as an active field of research and development. The underlying primary motivation is to judiciously exploit different flow topologies in a favorable way, either by focusing or attenuating the prevailing shock-wave systems. From a practical standpoint, design of a suitable barrier-configuration for efficient trapping of blast waves for shock-wave disaster prevention are often associated with complex flow geometries, consisting of multiparameter dependencies in order to achieve optimized shock mitigation. From the fundamental point of view, this problem involves basic structures of shock waves interacting with obstacles (such as cones, circular cylinders, spheres etc.), thus consisting of incident shock, direct or inverse Mach reflections, transmitted and reflected shocks, triple points and slip lines. This problem represents one of the simplest occurrences of strong shock/obstacle interactions and therefore an ideal test case for Computational Fluid Dynamics (CFD) solvers. In this respect, shock-wave diffraction phenomena over cones, cylinders and spheres had been studied by many researchers [1–7] in the past. Recent advances in computer sciences have led to increasing use of CFD methods as an important component of efficient process design with regard to complex fluid–solid interactions. Variety of numerical methods have evolved to deal with complex flow configurations. The immersed boundary methods (IBM) still remain popular as they can make use of relatively

⇑ Corresponding author. Tel.: +33 232 959 794; fax: +33 232 959 780. E-mail addresses: [email protected] (A. Chaudhuri), [email protected] (A. Hadjadj). 0021-9991/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2010.11.016

1732

A. Chaudhuri et al. / Journal of Computational Physics 230 (2011) 1731–1748

simple cartesian mesh configuration as a cost effective numerical tool. These methods allow to deal with complex geometries on rather simple cartesian meshes by introducing forcing terms on desired surfaces corresponding to the physical location of the complex boundaries. Numerous modifications and refinements have been proposed after Peskin [8] coined the term ‘‘immersed boundary method’’. Efficient implementation of boundary conditions is crucial while developing an immersed boundary algorithm. Detailed discussions have been presented by Mittal and Iaccarino [9] regarding the introduction of forcing functions to impose the boundary conditions through ‘‘continuous forcing approach’’ and ‘‘discrete forcing approach’’. It turns out that the direct-forcing approach bears less computational cost compared to other techniques [10]. On the other hand, hybrid cartesian immersed boundary approach (HCIB), the subsequent variant of IBM [11–16] with direct-forcing, has been evolved through a ghost-cell method and proved to be very promising for two- and three-dimensional computations. In their recent studies, Kang et al. [17,18] carried out direct numerical simulations using the immersed boundary approach. Although this method has been applied to a wide range of applications dealing with incompressible flow, studies with compressible flows [19,20] including strong shock waves are not abundant. In this paper, we report recent progress in the computation of the interaction of shock waves with rigid obstacles using high-order WENO scheme in conjunction with immersed boundary technique. The current research is motivated by the desire to develop reliable numerical tools to predict complex high-speed flows, especially for problems including shock waves interacting with instabilities, turbulence and solid obstacles. Subsequently, with more suitable numerical schemes and robust boundary methods, the proposed work will be able to shade more physical insight for the considered problems. The outline of this paper is as follows: Section 2 describes the governing equations as well as the numerical methods, including details of the flow solver and the IB technique. Section 3 describes the problem setup for supersonic flows past obstacles, including the description of the computed test-cases and the partial assessment of the proposed methodology through comparisons with analytical solutions from the shock wave theory. Further results of transient shock waves interacting with triangular prism (Schardin’s problem) and circular cylinder are presented and discussed. Finally, conclusions are drawn in Section 4.

2. Governing equations and numerical methods The Euler equations governing two-dimensional compressible inviscid flows can be written in conservative form as

@ t Q þ $F ðQÞ ¼ 0

ð1Þ

where Q is the conservative variable vector

Q ¼ ½q; q~ v ; qET

ð2Þ

and F ðQÞ indicates the inviscid fluxes

F ðQÞ ¼ ½q~ v ; q~ v ~ v þ PI; ðqE þ PÞ~ v T

ð3Þ

In the equations above: q; ~ v ; P and E are density, velocity vector, pressure and total energy per unit mass, respectively. The system is closed by the equation of state for ideal gas

P ¼ ðc  1ÞðqE  q~ v 2 =2Þ

ð4Þ

with the ratio of specific heats c = 1.4 for air. In the present work, the fifth-order WENO scheme [21], briefly described hereafter, is used to evaluate the numerical fluxes at cell interfaces. 2.1. WENO methodology The basic idea underlying the WENO approach consists of introducing a convex linear combination of low-order polynomial reconstructions that yields a high-order resolution in smooth regions and keeps the essentially non-oscillatory property near the discontinuities. For the sake of clarity, the WENO description is given for a one dimensional scalar problem. Exten1 sion to higher dimensions is straight forward. The WENO scheme computes numerical fluxes ^f  iþ1=2 at any cell interface, i þ 2, through reconstructed interpolated polynomials on r grid points per candidate stencil. The WENO candidate stencils can be expressed as Sk : {i  r + k + 1, i  r + k + 2, . . ., i + k}, for k = 0, . . ., r  1. To achieve high-order accuracy in smooth regions, the scheme uses a convex combination of these neighboring fluxes at i þ 12 and the interface approximation of the fluxes is given by

^f þ iþ1=2 ¼

r1 X

xk

k¼0

r1 X

þ arkl fðirþkþlþ1Þ

l¼0

where the nonlinear weights are

rk xk ¼ Pr

k¼0

rk

and

rk ¼

Xk ½ þ bk 2

ð5Þ

1733

A. Chaudhuri et al. / Journal of Computational Physics 230 (2011) 1731–1748

The smoothness indicators bk is given by k

b ¼

r1 r1 X X m¼0

!2 nrkml

þ firþkþlþ1

ð6Þ

l¼0

The values of the optimal coefficients arkl ; nrklm and Xk are reported in [21].  is a small parameter to prevent the denominator becoming zero and is usually taken as  = 106. For the evaluation of the numerical flux ^f  iþ1=2 , the procedure can be modified symmetrically with respect to xi+1/2. The total number of grid points available to the WENO scheme of Jiang and Shu [21] is 2r  1 with a maximum accuracy of 2r  1. 2.2. Immersed boundary method Regarding the IBM, the direct-forcing method is adopted in this study. The computational domain consists of fluid domain (Xfluid) and embedded solid domain (Xsolid) associated with the boundary function (s) wb (see Fig. 1). The ghost-cell finite difference approach is implemented to prescribe the desired boundary conditions on immersed boundaries. A typical ghost-fluid based method consists of first identifying the ghost points within the solid domain. Next, from each ghost point, a normal to the nearest boundary segment is built to locate the image of the corresponding ghost point into the fluid domain. The value of the flow-field variables at each image point is then interpolated from the surrounding fluid points. The flowfield variables at the image point are subsequently reflected back to the corresponding ghost point to ensure the requisite boundary conditions. The global cartesian grid can be represented as a set of ordered pair of coordinates as,

G ¼ fðxi ; yj Þ; i ¼ 1; . . . ; Nx ; j ¼ 1; . . . ; Ny g

ð7Þ

where Nx and Ny are the number of grid points in x and y directions, respectively. Depending on the geometry of the immersed solid(s), boundary function(s), wb, can be traced and the cartesian grid points belonging to fluid and solid can be identified. This information is stored as preprocessed data. Let us define a set of marker variables Zm as,

Zm ¼ ffi;j ; i ¼ 1; . . . ; Nx ; j ¼ 1; . . . ; Ny g where,

fi;j ¼



ð8Þ

1 if ðxi ; yj Þ 2 Xfluid 0

ð9Þ

otherwise

Therefore, the set of cartesian grid points belonging to fluid domain ðGfluid Þ and the set of cartesian grid points belonging to solid domain ðGsolid Þ can be given by,

Gfluid ¼ fðxi ; yj Þ 2 G s:t: fi;j ¼ 1g

ð10Þ

Gsolid ¼ fðxi ; yj Þ 2 G s:t: fi;j ¼ 0g

ð11Þ

Ωfluid NPk k

ψb

d

IP

j+1 BP j GP

j−1 y

Ω solid

x

i−1

i

i+1

Fig. 1. Immersed boundary treatment.

1734

A. Chaudhuri et al. / Journal of Computational Physics 230 (2011) 1731–1748

For Ns number of immersed solids, Gsolid ¼

G ¼ Gfluid

[

SNs

k¼1 Gsolid;k

so that the following relation holds

ð12Þ

Gsolid

Consequently, the solid points having at-least one fluid neighbor (ghost points GP in Fig. 1) can be identified so that the GP’s set can be given by,

GGP ¼ fðxi ; yj Þ 2 Gsolid if 9ðxk ; yl Þ 2 Gfluid for k ¼ i  2; . . . ; i þ 2; l ¼ j  2; . . . ; j þ 2g

ð13Þ

At this stage, one need to fix the flow-field parameters of the ghost points embedded inside the solid domain. In order to assign this, the image point (IP) is required to find for each ghost point orthogonal to the boundary, intersecting the boundary point (BP). Note that the definition of GGP includes at least two layers of ghost points, needed for the fifth-order WENO flux interpolation. This condition satisfies the requirement of conventional ghost points for the WENO scheme at computational boundaries. Each GP is thus associated with a unique IP, belonging to the fluid domain. The IP’s set can therefore be expressed as,

GIP ¼

n

o  h  i H H xH 2 Xfluid s:t: 9 an unique ðxi ; yj Þ 2 GGP and D xH ? wb ¼ D½ðxi ; yj Þ ? wb  i ; yj i ; yj

ð14Þ

where D½ðxa ; yb Þ ? wb  represents the orthogonal distance from a given point (xa, yb) to wb. The set of BP’s GBP can thus be defined as follows

^Þg GBP ¼ fð^x; y

ð15Þ

      H H ^ =2 for ðxi ; yj Þ 2 GGP and xH 2 GIP . For the sake of clarity, subscript ‘‘i, j’’ will be where ^ x ¼ xi þ xH i =2; y ¼ yj þ yj i ; yj omitted in the remaining text. The following linear approximation is used to impose either Dirichlet or Neumann boundary conditions for solid–fluid interface and to fix the flow-field variables at any GP,

QGP ¼ 2QBP  QIP

ð16Þ

In order to close the above equation, one need to fix the representative flow-field variables for each IP. This can be achieved by different interpolation techniques exploiting the available fluid domain flow-field variables of the four surrounding neighbor points NPk, k = 1, . . ., 4 (see Fig. 1). By applying a simple index searching algorithm to find a point entrapped inside four cell centers, each element of GIP can be associated with four NP’s. In this study, we have used the inverse distance weighting interpolation (IDW) [12,13] to evaluate QIP .

QIP ¼

4 X

ð17Þ

dk QNPk

k¼1

where,

dk ¼ gk

4 X

!1

gk

ð18Þ

k¼1 2

where gk ¼ 1=dk and dk’s are the distance of the corresponding IP from each NPk. Direct implementation of relation (18) may pose certain difficulties, particularly if NP k 2 Gsolid or dk ? 0 for a given k. These situations are illustrated in Fig. 2. For example, in Fig. 2(a) where NP 4 2 Gsolid , the contribution of NP4 should not be taken into account during interpolation. In order to exclude the contribution of neighbor points belonging to solid domain let us define ak in such a way that dk ? 1 if NP k 2 Gsolid as,

ak ¼



0 if NPk 2 Gsolid 1 otherwise

ð19Þ

Eq. (18) can thus be rewritten as,

dk ¼ ak gk

4 X

!1

ak gk

ð20Þ

k¼1

This simple numerical-fix preserves the convexity of the interpolation coefficients dk by automatically eliminating the contribution of NP k 2 Gsolid . For cases arising from dk ? 0 as shown in Fig. 2(b), contribution of NPl (l – k) can be neglected to enforce QIP ¼ QNPk (by Eq. (17)). In the present work, we assumed dk to be negligible if dk/d 6 1, with 1 = 1010 and pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d ¼ Dx2 þ Dy2 , where Dx, Dy are grid spacings. Flow-field data for ghost points can thus be pre-assigned before starting the space and time marching of the governing equations. It can be emphasized that the storage of geometrical data of ghost points are maintained separately for each space direction. Therefore, situation arising from sharp corner (see Fig. 2(c)) where any ghost point situated at a comparable distance from two space directions needs no special treatment. The flux in x-direction can be computed based on GP associated to IPx, whereas flux in y-direction can be computed utilizing GP associated with IPy. The convex interpolation thus preserves the unbiased stencil in both directions.

A. Chaudhuri et al. / Journal of Computational Physics 230 (2011) 1731–1748

Ωfluid

NP2

d2

d1

Ωfluid

NP1

IP d3

NPi

1735

IP

BP

NP3

NP4

BP

GP

Ωsolid

y x

Ωsolid

y x

GP

(a)

Ω fluid

(b) NP2

NP1 IPy BPy

NP3

IPx

NP4 BPx

y x

GP

Ωsolid (c)

Fig. 2. Schematic illustration of (a) one neighbor point belonging to solid domain, (b) image point very close to one neighbor point and (c) sharp edge with ghost points located at comparable distant in both directions.

Table 1 The WENO solver algorithm with immersed boundary method. STEP 1 STEP 2 STEP 3 STEP 4

STEP 5 STEP 6 STEP 7

? ? ? ?

? ? ?

generate grid data G create and store Zm using known boundary function(s) wb identify the set of ghost points GGP and the set of image points GIP based on Zm and wb (i) search neighbor points for each image points (ii) store the four neighbor points indices i, j (i = 1, . . ., Nx, j = 1, . . ., Ny) (iii) calculate and store dk, k = 1, . . ., 4 for each image point compute the numerical fluxes using the WENO solver update ghost point flow-field variables using Eqs. (16)–(20) before each RK steps of Eq. (21), based on requisite BC repeat from step 5 until the end of the simulation

Integrating Eq. (1) in time using an explicit third-order Runge–Kutta scheme together with the above direct-forcing scheme for immersed boundaries we get

8 H n n > < Q ¼ Q þ DtLðQ ÞZm QHH ¼ 14 ½3Qn þ QH þ DtLðQH ÞZm  > : nþ1 1 n ¼ 3 ½Q þ 2QHH þ 2DtLðQHH ÞZm  Q

ð21Þ

Since the solver developed in this study is fully explicit, one can assume stability of the overall integration under a Courant– Friedrichs–Lewy restriction, that is

Dt ¼ CFL=max½ðjuj þ cÞ=Dx þ ðjv j þ cÞ=Dy where c is the speed of sound, and u, v are velocity components. For all computations, the CFL number was fixed to 0.7. Finally, the main algorithm used for IBM implementation with the WENO solver is given in Table 1. In this study, the fifth-order WENO scheme is used for all reported simulations. 3. Results and discussion The proposed methodology is applied to compute two-dimensional quasi-steady and transient supersonic flows past obstacles in the limit of infinite Reynolds number. Extensive grid refinement tests were achieved to ensure that the global physical

1736

A. Chaudhuri et al. / Journal of Computational Physics 230 (2011) 1731–1748

features are well reproduced and the numerical results are mesh independent. Results presented thereafter contain approximately 1 (MP) Million of Points (for problems of Sections 3.1.1 and 3.1.2) and 8.4 MP (for problems of Sections 3.2.2 and 3.2.3), uniformly distributed with grid spacing of Dx = Dy  60 lm. Intensive parallel computing is involved in all simulations. Each simulation requires approximately 200–500 CPU hours using 24–40 processors on SGI ALTIX ICE machine. 3.1. Supersonic flow past solids The supersonic flow over a stationary solid obstacles is a typical problem that has been widely investigated both experimentally [22] and numerically [23–26]. Usually at high Reynolds number, the flow becomes unstable with large asymmetric recirculation zone and vortex shedding in downstream. The corresponding shedding frequency and the intensity of the vortex increase for higher Reynolds numbers. In this study, inviscid simulations are carried out for supersonic flow over triangular prism as well as circular cylinder at Mach 3.5 on a computational domain of lengths Lx  Ly = 120 mm  30 mm. Table 2 summarizes the simulation parameters for the two test cases. The supersonic inflow/outflow for inlet and outlet boundaries together with slip boundary conditions for top and bottom boundaries are fixed for these two test cases. Adiabatic walls with zero normal velocity (see Appendix A) conditions are enforced on solid boundaries. Initial conditions are considered as stagnant gas throughout the computational domain for these two cases. After an initial transient flow, a steady-state solution of supersonic stream is reached. In order to highlight the important effect of the interpolation functions presented in Section 2.2, the flow over triangular prism is simulated with interpolation functions (using Eqs. (16), (17) and (20)) and without interpolation functions: enforcing zero normal velocity based on IP which assumes the flow-field property of the nearest fluid point to be QIP ¼ QNPk if dk = min[di], i = 1, 4 for the kth neighbor. The obtained results are presented in the following subsections and compared to the analytical solutions of shock-wave theory. 3.1.1. Flow past triangular prism For flow over triangular prism, the deflection angle h is considered as 20° (see Fig. 3). From the following h–b–M relation

" tan h ¼ 2 cot b

2

M 2 sin b  1

# ð22Þ

M 2 ðc þ cos 2bÞ þ 2

the associated theoretical shock angle b is 34.58° for M = 3.5. Fig. 4 shows the magnified view of numerical Schlieren picture [27] for the solutions with and without interpolation functions. For simulation with interpolation, the calculated shock angle is found to be 34.5°, which is in excellent agreement with the theoretical value given above. Solution without interpolation and zero normal velocity condition at the solid surface yields not only reduced shock angle 31.9° (8% deviation from the theoretical value) but also generates instabilities and spurious oscillations near the boundary, which produces secondary shock waves. Solution with interpolation shown in Fig. 5, illustrates the whole flow field with multiple reflections of shock wave from the top and bottom boundaries, yielding complex shock/shock and shock/vortices interactions. Specifically, when the reflected shocks interact with the wake region, near the centerline of the flow, stream of weak vortices are generated and convected downstream. In order to analyze the flow characteristics and the vortex shedding in the wake region, the instantaneous pressure signal at xw = x/Lx  0.6 is recorded. Discussions on the analysis of pressure signal are presented after the description of the cylinder test-case.

Table 2 Geometrical and flow parameters for supersonic flow over solids, h: deflection angle, D: cylinder diameter and l: height of the triangular prism (dimensions are in mm). Case

M1

h (°)

D

l

Lx

Ly

Triangular prism Circular cylinder

3.5 3.5

20 –

– 5

11 –

120 120

30 30

Oblique shock M

β θ

Fig. 3. Schematic illustration of supersonic flow over a triangular wedge.

A. Chaudhuri et al. / Journal of Computational Physics 230 (2011) 1731–1748

1737

Fig. 4. Numerical Schlieren pictures for triangular prism with theoretical shock angle showed in red dash-dot line, (i) without interpolation and zero normal velocity at solid boundary (left) and (ii) with interpolation (right). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5. Numerical Schlieren picture of complete domain for triangular prism with interpolation.

3.1.2. Flow past circular cylinder Various numerical and experimental studies [28–32] showed that the supersonic flow around a circular cylinder, in particular the unsteady wake flow behind it contains many interesting features with a very complex structure. In this work, the challenge is to simulate this flow in great detail using IBM and high-order numerical schemes. The problem setup and the computational domain together with the boundary conditions are taken similar to the previous case. Unlike the triangular prism, a strong detached bow-shock is formed ahead of the cylinder. The local properties of the flow-field behind curved detached shocks allow the same type of analytic calculations as for the attached shock. From the numerical point of view, the two local angles h(s) and b(s) can be computed from the simulated flow-field data as well as the streamline properties along the shock curvature (s). Fig. 6 shows the bow shock and downstream unsteady structures in the whole domain. Computations are carried out for four different meshes to study the grid dependency, namely Mesh0: 504  126  0.061 MP, Mesh1: 1008  252  0.25 MP, Mesh2: 1424  354  0.5 MP, Mesh3: 2000  501  1 MP, Mesh4: 4000  1002  4 MP, and Mesh5: 8000  2004  16 MP. Fig. 7 shows the comparison of calculated h(s)–b(s)–M relation among the four meshes together with the relation from shock-wave theory. It can be seen that the progressive grid refinement reduces the deviation from the theoretical curve.

Fig. 6. Numerical Schlieren picture for supersonic flow past circular cylinder.

1738

A. Chaudhuri et al. / Journal of Computational Physics 230 (2011) 1731–1748

100

80

60

β 40 mesh4 mesh3 mesh2 mesh1

20

0

0

10

20

30

40

θ Fig. 7. Comparison of h–b–M plot for M = 3.5, —: shock wave theory.

For further comparison, the normalized distance between the shock and the  cylinder  leading edge are computed for different meshes and compared to the analytical formula given by Kx =D ¼ k1 exp k2 =M 21 , with k1 ’ 0.2 and k2 = 4.67 [33,34]. As shown in Table 3 when the mesh is refined, the numerical solution converges towards the theoretical value. Additionally, the surface pressure coefficient, C p ¼ 2ðP s  P1 Þ=q1 U 21 , is computed for different grids and results are shown in Fig. 8. It

Table 3 Comparison of bow-shock location for different meshes.

Kx/D

Theo.

Mesh5

Mesh4

Mesh3

Mesh2

Mesh1

Mesh0

0.293

0.30

0.30

0.30

0.31

0.33

0.34

Fig. 8. Surface pressure coefficient, inset: magnified view of flow domain containing the bow shock.

1739

A. Chaudhuri et al. / Journal of Computational Physics 230 (2011) 1731–1748

can be stated that grid independent solution is achieved for Mesh3 and Mesh4 for shock location and pressure surface properties. With regards to the convergence rate of the proposed method, it is worth mentioning that the presence of unsteady wake behind the obstacle, which includes acoustic waves and complex shock–vortex interactions, makes this problem quite difficult to statute on the formal accuracy of the scheme. However, in this paper we provided some measurements based on local and zonal L1 and L2 norms excluding the unsteady section, where surface angle is greater than 130 °. To assess the numerical error, the result of Mesh5 is used as a reference solution. Fig. 9 shows that for the present methodology the boundary accuracy is second order. Additionally, it is well known that the presence of strong shocks systematically degenerates the order of numerical schemes that have formally high-order accuracy in smooth region. The choice of WENO based schemes

Error

1

r=2

r=2 0.01

0.0001

0.125

Fig. 9. Local norms based on surface pressure coefficient (L1: N–N and L2: velocity ( ) and Dx0 represent the grid spacing of Mesh0.

10

1

Δx/Δx0

) and global norm L2 (pre-multiplied by factor 103) based on streamwise

6

Circular cylinder Triangular prism

-5/3

4

-1

Amplitude

10

10 8

2

P/Pinlet

10

6 4 2 0

10

0

10

0

0

0.1

0.2 0.3 time (ms) 10

1

0.4

0.5 2

10 Frequency (kHz)

10

3

Fig. 10. Comparison of power spectrum of the pressure signal at xw  0.6 (inset).

10

4

1740

A. Chaudhuri et al. / Journal of Computational Physics 230 (2011) 1731–1748

are advantageous due to their robustness for such situations. Keeping in view the application of the present methodology, we measured the zonal behavior of L2 norm to estimate the overall accuracy of the scheme. Fig. 9 depicts the global L2 norm measurement, showing that the present methodology preserves the accuracy higher than or equal to second order. The implementation of the present fifth-order WENO scheme without IBM has been also validated by solving isentropic vortex evolution problem (see Appendix B). It is evident that the maximum accuracy is achieved in smooth regions [35]. As mentioned before, the downstream pressure signals are recorded for the analysis of the flow evolution. Fig. 10 shows the pressure–time history for both triangular prism and circular cylinder cases. Although there exists a visual difference of these two raw pressure signal, it is useful to analyze these signals in frequency domain. The power spectrum of these Table 4 Geometrical and flow parameters for test cases with moving shocks, D: cylinder diameter, and the others parameters are defined in Figs. 4 and 11 (dimensions are in mm). Case

Ms

h (°)

D

b

Lx

Ly

Nx  Ny

1: Schardin’s problem 2: Circular cylinder

1.30 2.81

30° –

– 10

20 –

200 200

150 150

3360  2500 3360  2500

Lx

Ms

Ly

b

l y x

Incident shock wave Fig. 11. Schematic diagram of Schardin’s problem.

IS SL

TP1 MS

V

RS

TP2

Fig. 12. Different waves arising in Schardin’s problem with RS: reflected shock, IS: incident shock, SL: slip line, MS: Mach stem, TP: triple point and V: vortex.

1741

A. Chaudhuri et al. / Journal of Computational Physics 230 (2011) 1731–1748

pressure signals are shown in the same figure. It can be seen that there exits wide range of frequencies representing no specific dominant frequency of the signal. It can be recalled that the dissipation of mechanical energy is particularly dominating in the flow-field regions where the instantaneous velocity gradients are high. As a consequence, the energy dissipation is rather concentrated in the smallest turbulent eddies (close to the Kolmogorov length scale). Richardson visualized an inviscid cascade of energy from larger scales to smaller scales mainly driven by inertial forces. The main essence of Richardson’s cascade is that the energy transfer is indeed a multistage process. This involves a hierarchy of wide range vortices of varying size. It is obvious that viscosity plays no role in the energy cascade at high-Reynolds number, except at the smallest scales. Viscosity thus plays a passive role in Richardson’s cascade, mopping up whatever energy cascades down from above [36] and cannot influence the cascade itself. The spectrum for both cases shows less slope, scaling as f1 rather than f5/3. The

ray-shock theory (infinite wedge)

3

TP1

y/l

2

1 V

0

0

1

TP2

2

3

4

5

x/l Fig. 13. Comparison of triple point trajectory and locus of vortex center ( : experimental data [7],

: numerical results [7], – : present simulation).

Fig. 14. Comparison of numerical shadowgraph (top-half) and experimental shadowgraph [7] (bottom-half), view dimensions: 50–137 mm in x-direction and 75–106 mm in y-direction.

1742

A. Chaudhuri et al. / Journal of Computational Physics 230 (2011) 1731–1748

absence of viscous terms in the present computation is manifested by the presence of high amplitude, associated with high frequencies of the power spectrums. It is also worth mentioning that the existence of the interaction between transverse shocks with the convected vortices in downstream influences the cascade law as it affects the vortex merging process, commonly observed in usual 2D turbulence. 3.2. Moving-shock/obstacle interactions In this section, we consider the case of non-stationary shock waves past obstacles. This kind of situation arises, for example, when a valve or membrane is suddenly opened and the resulting shock hits an obstacle placed in its passage while propagating downstream. With the advent of the shock front, an unsteady flow is induced around the obstacle accompanied by complex wave patterns, produced by multiple reflections and interactions. In the course of these processes, flow separation (for bluff obstacles) and/or chocking (for large-blockage ones) are considered to occur. Generally, the transmitted shock is

Fig. 15. Evenly spaced snapshots of shock–vortex dynamics for Schardin’s problem, AW: accelerated wave, DW: decelerated wave (see Fig. 12 for other notations), view dimensions: 47–130 mm in x-direction and 41–116 mm in y-direction.

1743

A. Chaudhuri et al. / Journal of Computational Physics 230 (2011) 1731–1748

weaker than the incident one and is retarded with respect to the latter. The flow-field suffers momentum loss due to the resistance of the obstacle. To investigate such complex phenomena, two simple-shaped obstacles are considered, namely the triangular prism and circular cylinder cases. Comparisons are made with previous experimental and numerical results of literature [1,3–5,7,6]. Table 4 summarizes the detail of the simulation parameters for the two problems. The computational domain is chosen sufficiently large to avoid any inflow/outflow wave reflections or boundary-layer interaction other than the solid obstacle. 3.2.1. Flow initialization In order to initialize the flow-field with a moving shock, Rankine–Hugoniot relations are applied [37]. For a normal shock propagating at a constant speed, the velocity of the moving fluid is given by

uL 1 ¼ ðn  1Þ=M s cR c s

ð23Þ

where Ms is the shock Mach number, c isthe ratio  of specific heats and ns is the pressure jump across the shock, also known as shock strength, defined as ns ¼ M2s þ j M 2s  1 , with j = (c  1)/(c + 1). Subscript ‘‘L’’ denotes the left shocked state of the moving flow, while subscript ‘‘R’’ characterizes the stagnant downstream flow-field. The relationships between the pressure jump across the shock and the density and temperature ratio, respectively are given by

qL =qR ¼ ð1 þ jns Þ=ðj þ ns Þ T L =T R ¼ ns ðj þ ns Þ=ð1 þ jns Þ

ð24Þ ð25Þ

From known values of Ms, PR, TR using the above set of equations, the flow-field variables are initialized for the two simulations presented hereafter. 3.2.2. Shock/triangular prism interaction Fig. 11 illustrates the flow configuration of Schardin’s problem [1]. The static pressure and temperature of the lowpressure side are taken 0.05 MPa and 300 K, respectively. A schematic diagram of various waves appearing in Schardin’s problem is depicted in Fig. 12. The flow characteristics consist of unsteady evolution of multiple Mach stems, triple points, reflected (attenuated and accelerated) shocks, slip-lines and vortices. The reader is referred to [7] for further details and description of this test-case. Fig. 13 illustrates the comparison of present simulation and results by Chang and Chang [7] in terms of triple point trajectory and vortex location. It can be seen that the simulated time evolution of the two triple points (TP1 and TP2 ) and the locus of vortex core (V) are comparable to the experimental and numerical results of Chang and Chang [7]. Instantaneous numerical shadowgraph picture at t = 138 ls is compared with the experimental shadowgraph reported in [7]. It is evident from Fig. 14 that the essential flow characteristics are reproduced by the simulation with high fidelity.

4

3

y/D

TP1 2

1

TP2 0

0

1

2

3

4

5

6

x/D Fig. 16. Comparison of triple point trajectory: case 2 ( : experimental data [4],

: experimental data [3],

: numerical results [6], – : present simulation).

1744

A. Chaudhuri et al. / Journal of Computational Physics 230 (2011) 1731–1748

Fig. 17. Numerical Schlieren picture (left) and streamlines (right) at t = 30 ls (see Fig. 12 for notations), view dimensions: 38–100 mm in x-direction and 48–102 mm in y-direction.

Fig. 15 shows a sequence of snaps of the flow evolution. It can be clearly seen that the inviscid computations are sufficient enough to capture the global flow characteristics and results are in excellent agreement with the experimental data.

3.2.3. Shock/circular cylinder interaction The last example concerns the shock/cylinder interaction problem which is also one of the well studied test cases in the literature [2,3,5,6]. The flow consists of a planar shock moving at M = 2.81 and impinging on a solid cylinder. The computational domain is taken similar to the Schardin’s case (Table 2). The computed trajectories of the two triple points are depicted in Fig. 16. The comparisons appear substantially satisfactory for both upper and lower shock wave systems. The detail of complex wave interactions is further illustrated in Fig. 17. The flow pattern appeared symmetric with respect to the centerline which is similar to that observed in Schardin’s problem. The evolution of diffraction process including regular reflection, transition to Mach reflection, complex wake/shock and shock/shock interactions are correctly reproduced by the present simulation. The formation of a pair of counter-rotating vortices together with the appearance of a curved reflected shock can be seen from the streamlines. The front and rear stagnation points as well as flow separation are also clearly evident from Fig. 17.

4. Conclusions This paper provides a numerical methodology for the solution of two-dimensional shock/obstacle interaction problems using a high-order WENO scheme in conjunction with an immersed boundary method. The implementation of IBM is based on direct-forcing and ghost-fluid approach. This method is particularly suited for computing complex flows of practical interest. An interpolation algorithm is developed for more stable boundary conditions with simple implementation techniques. Detailed discussions of the method are presented on the interpolation of flow variables, direct-forcing of ghost cells, resolution of immersed-boundary points and internal treatment. A simple numerical-fix is proposed to preserve the convexity of the weighting coefficients for interpolation functions. The developed numerical solver is validated against well documented steady and unsteady supersonic test problems with 5th order classical WENO method and IBM. It is evident that inappropriate estimation of flow variables for embedded ghost points leads to erroneous numerical predictions. The computations without interpolation functions for the supersonic flow over triangular prism is carried out to demonstrate the associated inaccuracy incurred. The implementation of the interpolation functions together with the proposed numerical-fix comprehensively validated against the solution of oblique shock wave theory. Error analysis based on the supersonic flow over cylinder reveals that the local accuracy of the proposed IBM strategy is second order, whereas the global L2 norm indicates an accuracy of higher than or equal to second order. Keeping in view that the accuracy of high-order schemes degenerates in presence of shock, the choice of high-order WENO schemes undoubtedly favorable due to its robustness. As well as the behavioral study of modified WENO schemes in presence of shock/turbulence interaction is yet an active field of research. In this aspect, recent study [38] reveals the advantage of new strategies with balanced filter based modification of WENO schemes. In this conjecture, it is worth to mention that the present study can be further extended using various modified WENO schemes coupled with the proposed IBM to explore its effect on global accuracy and robustness. Additionally, two different test cases with moving shocks, bearing complex shock–vortex dynamics are used to assess the performance and the robustness of the method for transient supersonic flows. The comparison with the experimental results clearly shows once more the performance of the proposed algorithm. Although the current results agree well with the analytical solutions as well as the experimental data, inclusion of viscous terms are indeed essential for more realistic predictions of complex flows.

A. Chaudhuri et al. / Journal of Computational Physics 230 (2011) 1731–1748

1745

Acknowledgments This work was performed using HPC resources from GENCI [CCRT/CINES/IDRIS] (Grant 2009-0211640). The authors are very much thankful to Dr. Se-Myong Chang from Kunsan National University (Korea) for providing the experimental shadowgraph for comparison and validation. Appendix A. Implementation of zero-normal velocity at the solid body Let us consider any arbitrary surface as shown in Fig. A.18. The following steps can be followed to implement the zero normal velocity at the solid boundary and to fix the velocity components at the associated GP. Suppose at any IP, vxIP, vyIP are the velocity components in the axis directions. On the other hand, vn and vs are the components along surface normal and tangent directions. 1. Calculate hs, the angle between x-axis and velocity vector ~ v of any IP

hs ¼ tan1



v yIP v xIP

ðA:1Þ

2. Based on the quadrant in which hs belongs, modify hs as follows

hs :¼ ðm þ nÞP þ hs

ðA:2Þ

where

( m¼

0

if 0 6 hs 6 12 p

1 otherwise

( ;

n¼

1 if 0

3 2

p 6 hs 6 2p

otherwise

3. Calculate as, the angle between x-axis and the surface normal direction gs as follows

as ¼ tan1



yIP  yGP x IP  xGP

ðA:3Þ

4. Modify as similar to hs in step (3) 5. Calculate cs, vxGP and vyGP as follows

cs ¼ P þ 2as  hs ; v xGP ¼ j~ vj cos cs ; v yGP ¼ j~ v j sin cs

ðA:4Þ

Appendix B. Error analysis for WENO scheme without IBM We consider the evolution of a 2D inviscid isentropic vortex as a test problem without shocks or turbulence to assess the performance of the WENO scheme used in the present work without IBM. Several authors [39–41] presented the similar test

Fig. A.18. Schematic view for implementation of zero normal gradient of velocity at the solid surface.

1746

A. Chaudhuri et al. / Journal of Computational Physics 230 (2011) 1731–1748

case for testing numerical schemes with respect to the vortex preservation. This test case consists of a uniform twodimensional flow field, in which an isentropic vortex is created by the following perturbation

C

ðdu; dv Þ ¼

2 Þ=2

2p

eð1^r

; xÞ; ðy

dT ¼ 

ðc  1ÞC2 1^r2 e 8cp2

Þ ¼ ðx  x0 ; y  y0 Þ. x0 and y0 are the initial center of the vortex core, and where C = 5 is the vortex strength and ð x; y ^r2 ¼  2 . Since the flow is isentropic and the gas is assumed perfect p and q satisfies p/qc = 1. The initial conditions x2 þ y are taken as follows

v ¼ v 1 þ dv ;

u ¼ u1 þ du;

T ¼ T 1 þ dT;

1

q ¼ T c1

where p1 = 1, T1 = 1, u1 = v1 = 0. Here, the perturbation induced by the vortex is not large enough to produce a significantly nonlinear effect. The exact solution with a given initial states is just a passive convection of the vortex with the mean velocity. Thus, the solution should remain unchanged as time evolves. Initially, a vortex with a radius ^r ¼ 5 is placed at the center of the computational domain covering (x, y) = (10, 10). The boundary conditions are taken periodic in both directions. The CFL number is set to 0.8 and the results are shown at normalized time equal to 1. We measured the relative and absolute errors during time evolution to estimate the order of convergence. If the method is of rth order, the order of convergence, r, for a uniform mesh is then given by

h r ¼ ln

E Nn =E Nn

0

i ðB:1Þ

lnðN0 =NÞ 0

where E Nn and E Nn are the errors in Ln norms using the mesh of grid points N  N and its double N 0  N 0 , respectively. At t = 1, the computed vortex is compared with the exact solution and the convergence rate are calculated according to Eq. (B.1).

1

0.9

0.8 N = 20 N = 40 N = 80 N = 160 N = 320 N = 640 N = 1280

ρ 0.7

0.6

0.5

3

4

5

6

7

x Fig. B.19. Grid refinement analysis. Density profiles at t = 1.

Table B.5 Convergence rates of WENO5. N

Relative L1 error

r

L1 error

r

20 40 80 160 320 640 1280

4.9E2 5.2E3 2.0E4 4.2E6 1.2E7 3.2E9 8.6E11

– 3.2 4.7 5.6 5.1 5.2 5.2

1.4E1 1.9E2 9.3E4 1.6E5 5.0E7 1.3E8 1.0E9

– 2.9 4.3 5.9 5.0 5.3 3.7

A. Chaudhuri et al. / Journal of Computational Physics 230 (2011) 1731–1748

1747

Fig. B.20. Relative L1 ( ) and L1 ( ) of WENO scheme. Inset: Density contour for most refined mesh.

Fig. B.19 shows the influence of the grid resolution on the accuracy of the solution based on density profile obtained for different meshes. Table B.5 summarizes the measure of errors compared to the exact solution. It is clear from Table B.5 and Fig. B.20 that the WENO scheme used in the present work preserves a fifth-order accuracy. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30]

H. Schardin, High frequency cinematography in the shock tube, J. Photo Sci. 5 (1957) 19–26. G.B. Whitham, A new approach to problems of shock dynamics, Part I: Two-dimensional problems, J. Fluid Mech. 2 (1957) 145–171. A.E. Bryson, R.W.F. Gross, Diffraction of strong shocks by cones, cylinders and spheres, J. Fluid Mech. 10 (1961) 1–16. J. Kaca, An interferometric investigation of the diffraction of a planar shock wave over a semicircular cylinder, UTIAS Technical Note, 269, 1988. J.Y. Yang, Y. Liu, H. Lomax, Computation of shock wave reflection by circular cylinders, AIAA J. 25 (5) (1987) 683–689. J. Zoltak, D. Drikakis, Hybrid upwind methods fir the simulation of unsteady shock-wave diffraction over a cylinder, Comput. Methods Appl. Mech. Eng. 162 (1998) 165–185. S. Chang, K. Chang, On the shock–vortex interaction in Schardin’s problem, Shock Waves 10 (2000) 333–343. C.S. Peskin, Flow patterns around heart valves: a digital computer method for solving the equations of motion, PhD Thesis, Physiol., 378, Albert Einstein Coll. Med., Univ. Microfilms, 1972, pp. 72–30. R. Mittal, G. Iaccarino, Immersed boundary methods, Annu. Rev. Fluid Mech. 37 (2005) 239–261. G. Iaccarino, R. Verzicco, Immersed boundary technique for turbulent flow simulations, Appl. Mech. Rev. 56 (2003) 331–347. E.A. Fadlun, R. Verzicco, P. Orlandi, J. Mohd-Yusof, Combined immersed finite-difference methods for three-dimensional complex flow simulations, J. Comput. Phys. 161 (2000) 35–60. Y. Tseng, J.H. Ferziger, A ghost-cell immersed boundary method for flow in complex geometry, J. Comput. Phys. 192 (2003) 593–623. T. Gao, Y. Tseng, X. Lu, An improved hybrid cartesian/immersed boundary method for fluid–solid flows, Int. J. Numer. Methods Fluids 55 (2007) 1189– 1211. A. Dadone, B. Grossman, Ghost-cell method for inviscid two-dimensional flows on cartesian grids, AIAA J. 42 (12) (2004) 2499–2507. A. Gilmanov, F. Sotiropoulos, A hybrid cartesian/immersed boundary method for simulating flows with 3D, geometrically complex moving bodies, J. Comput. Phys. 207 (2005) 457–492. A.F. Shinn, M.A. Goodwin, S.P. Vanka, Immersed boundary computations of shear- and buoyancy-driven flows in complex enclosures, Int. J. Heat Mass Transfer 52 (2009) 4082–4089. S. Kang, G. Iaccarino, F. Ham, DNS of buoyancy-dominated turbulent flows on a bluff body using the immersed boundary method, J. Comput. Phys. 228 (2009) 3189–3208. S. Kang, G. Iaccarino, F. Ham, P. Moin, Prediction of wall-pressure fluctuation in turbulent flows with an immersed boundary method, J. Comput. Phys. 228 (2009) 6753–6772. R. Ghias, R. Mittal, H. Dong, A sharp immersed boundary method for compressible viscous flows, J. Comput. Phys. 225 (2007) 528–553. Q. Liu, O.V. Vasilyev, A Brinkman penalization method for compressible flows in complex geometries, J. Comput. Phys. 227 (2007) 946–966. G. Jiang, C.W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1996) 202–228. C.-S. Kim, Experimental studies of supersonic flow past a circular cylinder, Journal of the Physical Society of Japan 11 (1956) 439–445. P.W. Duck, The inviscid axisymmetric stability of the supersonic flow along a circular cylinder, J. Fluid Mech. 214 (1990) 611–637. A. Kluwick, P. Gittler, R.J. Bodonyi, Viscous–inviscid interactions on axisymmetric bodies of revolution in supersonic flow, J. Fluid Mech. 140 (1984) 281–301. M. Pandolfi, F. Larocca, Transonic flow about a circular cylinder, Comput. Fluids 17 (1989) 205–220. P. De Palma, M.D. de Tullio, G. Pascazio, M. Napolitano, An immersed-boundary method for compressible viscous flows, Comput. Fluids 35 (2006) 693– 702. A. Hadjadj, A. Kudryavtsev, Computation and flow visualization in high speed aerodynamics, J. Turbul. 6 (16) (2005) 33–81. Y.A. Panov, A.I. Shvets, Pressure oscillation on a plane upstream of an obstacle in a supersonic flow, Fluid Dyn. 33 (1998) 56–60. V.S. Murthy, W.C. Rose, Form drag, skin friction, and vortex shedding frequencies for subsonic and transonic gross flow on circular cylinder, AIAA Paper, No. 77687, 1977. J.M. Macha, Drag of circular cylinders at transonic Mach numbers, J. Aircraft 14 (1977) 605.

1748

A. Chaudhuri et al. / Journal of Computational Physics 230 (2011) 1731–1748

[31] V.A. Bashkin, A.V. Vaganov, I.V. Egorov, D.V. Ivanov, G.A. Ignatova, Comparison of calculated and experimental data on supersonic flow past a circular cylinder, Fluid Dyn. 37 (2002) 473483. [32] M.M. Zdravkovich, Flow Around Circular Cylinders, Applications, vol. 2, Oxyford University Press, 2003. [33] F.S. Billig, Shock-wave shapes around spherical and cylinder nosed bodies, J. Spacecraft Rockets 4 (1967) 822–823. [34] O. Boiron, G. Chiavassa, R. Donat, A high-resolution penalization method for large Mach number flows in the presence of obstacles, Comput. Fluids 38 (2009) 703–714. [35] S. Gottlieb, D. Gottlieb, C.W. Shu, Recovering high-order accuracy in WENO computations of steady-state hypersonic systems, J. Sci. Comput. 28 (213) (2006) 307–318. [36] P.A. Davidson, Turbulence An Introduction for Scientists and Engineers, Oxford University Press, 2004. [37] J.D. Anderson, Modern Compressible Flow with Historical Perspective, McGraw-Hill Book Company, 1982. [38] W. Wang, H.C. Yee, B. Sjögreen, T. Magin, C.W. Shu, Construction of low dissipative high-order well-balanced filter schemes for non-equilibrium flows, J. Comput. Phys. (2010), doi:10.1016/j.jcp.2010.04.033. [39] D. Gottlieb, E. Turkel, Dissipative two-four methods for time dependent problems, Math. Comput. 30 (1976) 703. [40] F. Davoudzadeh, H. McDonald, B.E. Thompson, Accuracy evaluation of unsteady CFD numerical schemes by vortex preservation, Comput. Fluids 24 (1995) 883. [41] H.C. Yee, N.D. Sandham, M.J. Djomehri, Low-dissipative high-order shock-capturing methods using characteristic-based filters, J. Comput. Phys. 150 (1999) 199–238.