Characteristics of overexpanded nozzle flows in imposed oscillating condition

International Journal of Heat and Fluid Flow 46 (2014) 70–83

Contents lists available at ScienceDirect

International Journal of Heat and Fluid Flow journal homepage: www.elsevier.com/locate/ijhff

Characteristics of overexpanded nozzle flows in imposed oscillating condition A.B.M. Toufique Hasan ⇑ Department of Mechanical Engineering, Bangladesh University of Engineering and Technology (BUET), Dhaka 1000, Bangladesh

a r t i c l e

i n f o

Article history: Received 23 October 2013 Received in revised form 6 January 2014 Accepted 6 January 2014 Available online 31 January 2014 Keywords: Supersonic nozzle Overexpanded flow Shock waves Mach stem Shock separation

a b s t r a c t A computation study is performed to investigate the effect of imposed oscillation of nozzle pressure ratio (NPR) on the flow structure in a two-dimensional, non-axisymmetric supersonic converging–diverging nozzle. In this study, the overexpanded flow conditions are considered which are dominated by the shock-induced boundary-layer interaction and corresponding free shock separation. The computational results are well validated with the available experimental measurements. Results showed that the internal flow structure of the nozzle is dependent on the process of change of pressure ratio during the oscillation. Distinct flow structures are observed during increasing and decreasing processes of the change of pressure ratio even when the nozzle is at the same NPR. Irreversible behaviors in the locations of free shock separation, Mach stem, and the strength of Mach stem are observed at the same NPRs during this oscillation. However, the nozzle thrust performance does not show the same order of irreversibility as in the cases of shock structures. Further, the effect of oscillation frequency is explored on this irreversible behavior. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction The research on flow through a supersonic converging–diverging nozzle is of great importance in many advanced engineering applications, especially in the aeronautical and aerospace industries. When a supersonic nozzle operates at pressure ratios below its design point (over-expanded), shock waves form inside the nozzle. In these cases, shock waves interact with the boundary layer and as a consequence flow separation occurs (Gatski et al., 2013). However, this separation often grouped into two primary categories namely-free shock separation (FSS) and restricted shock separation (RSS) (Frey and Hagemann, 1999; Hadjadj and Onofri, 2009). In FSS, the separation region downstream of the shock fails to reattach for the remaining length of the nozzle. However, the RSS is characterized by a small separation region or ‘‘bubble’’ which exists immediately downstream of the shock wave. In this region, the mean flow circulates and separates or tilts away from the wall before the flow reattaches and continues down the length of the nozzle as an attached boundary layer. The compressible flow systems can undergo pulsating conditions which are observed in many real field applications such as in airblast nozzle (injector), Rijke tube burners (McIntosh and Rylands, 1996) and so on. Strasse (2011) investigated the characteristics of the three-stream coaxial airblast nozzle which is used ⇑ Tel.: +880 173 071 4444. E-mail address: toufi[email protected] 0142-727X/$ - see front matter Ó 2014 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.ijheatfluidflow.2014.01.001

to generate an atomized fuel stream for a large-scale reactor. A pulsatile spray pattern was captured in the experiments and computations with frequencies between 75 Hz and 600 Hz. Depending on the frequencies, three distinct flow patterns was observed namely- spray sheds at outer air/film interface, spray sheds at the inner air/film interface and the bulk flapping of annular liquid film. Further, similar pulsating condition is observed in the interaction of compressible jets in supersonic cross-flow (Semlitsch et al., 2013). In the inviscid computation of Semlitsch et al. (2013), it was documented that the shock pattern is greatly influenced by the jets and flow field becomes pulsating with dominating frequency of 266 Hz (Strouhal-number of 0.5). The complex unsteady shock behavior in a simple over-expanded planar nozzle has been studied by Johnson and Papamoschou (2010). The experiments have revealed a large asymmetry in the two-dimensional shock structure as well as in the region of separation downstream of the shock. Large scale unsteadiness of the shock wave was also present. To better understand the flow physics of asymmetric nozzle flows, Xiao et al. (2009) performed an unsteady Reynolds-averaged Navier–Stokes calculation on the same geometry. Agreement in the mean shock location, pressure profile, and exit velocity profile was surprisingly good in 2D solutions. Several turbulence models were tested which produced slightly different results and all were able to capture the shock wave asymmetry. Further, in the experiments of Papamoschou et al. (2009), the shock was in lambda structure on one side wall and unsteady. However, there was no evidence of resonant tones.

A.B.M. Toufique Hasan / International Journal of Heat and Fluid Flow 46 (2014) 70–83

71

Nomenclature A a0 E F FSS f ht k M _ m NPR p p01 pb R RSS SBLI T t T01 u

area, m2 speed of sound, m/s total energy, kJ nozzle thrust, kN free shock separation frequency of NPR oscillation, Hz half-height at nozzle throat, m turbulent kinetic energy, m2/s2 Mach number mass flow rate, kg/s nozzle pressure ratio static pressure, kPa inlet total pressure, kPa back pressure, kPa universal gas constant, J/kg K restricted shock separation shock-boundary layer interaction time period of NPR oscillation, ms/temperature (K) time, s inlet total temperature, K x-velocity component, m/s

Time-resolved wall pressure measurements indicated that the shock oscillates in a piston-like manner and most of the energy of the oscillations is at low frequency. Recently, a mechanism for unsteady separation and shock motion in over-expanded planar nozzles has been identified by Olson and Lele (2013) using numerical simulation and a reduced order model. A series of calculations have been conducted to elucidate the effects of Reynolds number, nozzle geometry, and shock strength on this unsteadiness. It was found that the higher Mach number shock waves produce larger oscillations. However, the incoming boundary layers that are at higher Reynolds numbers produce oscillations with smaller amplitudes. The frequency is primarily dependent on the distance from the mean shock location to the exit of the nozzle and the shock Mach number. Further, the results suggested that the mechanism for the unsteadiness is a feedback loop between the separated shear layer and the shock wave. Beside these, the response of a transonic channel flow when the shock-wave is subjected to a periodic motion at a well defined frequency was investigated by Bur et al. (2006). Two component laser Doppler velocimeter probing was carried out for a shock oscillation frequency of 40 Hz. Phase-averaged field results enabled both the wave propagation in the core flow and the response of the boundary layer subjected to an oscillating shock wave. It was concluded that in the shock oscillation region, no phase lag was observed between velocities in the core flow and in the boundary layers whereas a significant one has appeared in the downstream subsonic region. An experimental study of force oscillation on normal shock wave in a parallel-walled duct has been conducted by Bruce and Babinsky (2008). The frequencies were varied from 16 Hz to 90 Hz. Normal shocks have been observed to undergo oscillatory motion in response to an imposed varying pressure ratio. It concludes that the mechanism by which shocks respond to back pressure variations was due to change their relative strength by moving so that their relative Mach number matches the pressure jump. These changes in relative shock strength can lead to changes in the extent of boundary layer separation and SBLI structure. However, in recent years, hysteresis phenomena in fluid flow systems drew attention for their great variety of industrial and engineering applications. Ben-Dor et al. (2003) performed a numerical simulation to investigation the flow-Mach-number-induced

v x y

y-velocity component, m/s location in x-direction, m location in y-direction, m

Greek symbols e turbulence dissipation rate, m2/s3 q density, kg/m3 c specific heat ratio l viscosity, Pa s s viscous stress tensor x specific dissipation rate, s1 Subscripts e exit eff effective i ideal IS incident shock wave MS Mach stem s shock wave t throat/turbulent

hysteresis phenomenon in the shock-on-shock interaction of conical shock waves. The specific geometry of the curvilinear cone was chosen in order to promote the RR M MR (Regular Reflection M Mach Reflection) transition. The upstream Mach number is varied from 2.5 to 5.0. As a results, a multiplicity of hysteresis loops was observed. It was concluded that there are flow Mach number ranges in which four different wave configurations, three inviscid and one viscous can be obtained for identical flow conditions. The different wave configurations for identical flow Mach numbers were associated with different pressure distributions. This study further complements an earlier study by Ben-Dor et al. (2001) in which an angle of-incidence-induced hysteresis was investigated. In the experiments, the angle of interaction between the two conical shock waves was altered continuously by moving the curvilinear cone along the axis of symmetry and as a result both the RR ? MR and the MR ? RR transitions were observed and recorded. This viscous dependent hysteresis loop was also obtained numerically by Burstchell et al. (2001) in their Navier–Stokes simulations. In case of planar overexpanded jets, the similar hysteresis behavior was observed by Hadjadj et al. (2004). The dependencies of the flow pattern, and the pressure distribution, inside an air intake on the variations of flight speed of a supersonic/hypersonic aircraft should be accounted for in the designing of intakes and flight conditions for supersonic and hypersonic vehicles. Onofri and Nasuti (2001) first related the hysteresis phenomenon in the RR M MR transition to the performance of air intakes. However, their numerical study was performed using a two-dimensional geometry, which resembles actual supersonic air intakes. Moreover, Shimshi et al. (2011) performed a viscous simulation to investigate the flow in overexpanded nozzles of planar and tapered shapes at different nozzle pressure ratio (NPR). The simulation reveals complex changes in the flow structure as the ratio between the ambient and the stagnation pressures is increased and decreased. Further, the simulation in planar nozzles showed the existence of a hysteresis process in the transition from regular to Mach and from Mach to regular reflections inside an ideal nozzle. However, such a process did not appear in tapered nozzles. Moreover, the transient side-load phenomena in a rocket nozzle during start-up process was investigated by Lijo et al. (2010). In this case, a hysteresis phenomena of shock structure from FSS to RSS was observed.

72

A.B.M. Toufique Hasan / International Journal of Heat and Fluid Flow 46 (2014) 70–83

In the present study, a numerical computation is performed to investigate the flow structure inside an overexpanded nozzle in response to oscillating flow conditions. The NPR is oscillated sinusoidally with time. Flow structures are captured in increasing and decreasing processes during NPR oscillation. Moreover, the effect of oscillation frequency is also studies on the flow structures. Computed results of separated nozzle flows are first verified with the available experimental results at steady NPRs. 2. Governing equations and computational methods 2.1. Governing equations The flow through a converging–diverging (C–D) nozzle is governed by time dependent, compressible, viscous form of the fluid flow conservation equations. For simulating the compressible flow, the total energy equation including viscous dissipation is also included and coupled to the set with the perfect gas law. The thermodynamics and transport properties for air are held constant; their influence was not found to be significant. The governing equations can therefore be written in the Cartesian indices notation as follows: Continuity:

@ q @ðqui Þ þ ¼0 @xi @t

ð1Þ

Momentum:

@ðqui Þ @ðqui uj Þ @p @ sij þ ¼ þ @t @xj @xi @xj

ð2Þ

Energy:

  @ðqEÞ @ðui ðqE þ pÞ ~ @T þ ¼ r  aeff þ uj ðsij Þ @t @xi @xi

ð3Þ

Gas-law:

q¼

p RT

ð4Þ

with

sij ¼ leff

  @ui @uj 2 @uk  leff þ dij 3 @xj @xi @xk

ð5Þ

In the above equations E is the total energy which is equal to h  p/

q + (u2 + v2)/2 and h is the sensible enthalpy. leff is the effective viscosity in the flow field and it is the sum of absolute viscosity (laminar), ll and turbulent viscosity, lt. To compute the homogenous turbulent contributions to momentum and energy transport in the present flow field, a two-equation k–x shear stress transport (SST) turbulence model of Menter (1994) is included with the above equations. This model uses a smooth blending between the standard k-e model in freestream and the k–x model of Wilcox (1986) near the wall. Further, this model formulates the modified turbulent viscosity to account for the transport effects of the principal turbulent shear stress. This model was used in various researches on compressible flow systems and captured the turbulence with good accuracy (Balabel et al., 2011; Strasse, 2011). 2.2. Computational methods The governing equations as mentioned in Section 2.1, are first transformed from physical space to computational space by generalized coordinate transformation to enhance the efficiency and accuracy of numerical scheme and to simplify the implementation of boundary conditions. The resultant Jacobian transformation of

the coupled system of equations is given in Appendix A. These equations are discretized using a control volume technique. For all equations, convective terms are discretized using a second-order upwind scheme, inviscid fluxes are derived using a second order flux splitting (Roe, 1981, 1986), achieving the necessary upwinding and dissipation close to shock waves. The interface flux is determined by separate terms, which depends on the upstream and downstream sides of the face, so that the information passed through the face contains the flow characteristics. Diffusion terms are always cast into a central difference form. To reconstruct the gradients of variables, a least square cell based method is used in the present numerical simulation. A TVD slope limiter (Barth and Jespersen, 1989) is implemented which uses the Minmod function to limit and restrict the reconstructed solution overshoots and undershoots on the cell faces. The resulting system is time-preconditioned, in order to overcome the numerical stiffness encountered at low Mach number. The discretized system is solved in a coupled way with a Block Gauss–Seidel algorithm. For the time derivatives, an implicit multistage time stepping scheme, which is advanced from time t to time t + Dt with a second order Euler backward scheme for physical time and implicit pseudo-time marching scheme for inner iteration, is used. The criterion for assessing convergence was based on the root mean square of the resides. The computations are stopped when the residual of mass, momentum, energy, and turbulence properties fall below 1  105. Since a time marching solver is used, the required time step is set owing to a Courant–Friedrichs–Lewy (CFL) condition that ensures the acoustic waves and flow physics are properly tracked by the solver. Acoustic velocity scale is used to determine the CFL criteria. During the iterations, the CFL is set to be 0.9. 2.3. Computational domain and boundary conditions In the present study, the geometric dimension of the supersonic nozzle is taken from Hunter (2004) as shown in Fig. 1a. The test nozzle is a subscale, non-axisymmetric, two-dimensional C–D nozzle with an expansion ratio Ae/At = 1.797. Based on one-dimensional theory, the nozzle has a design nozzle pressure ratio (NPR) of 8.78 and an exit Mach number of Me = 2.07. The origin of (x, y) coordinates is located at inlet of the nozzle. The half height of nozzle throat, ht is 27.48 mm. The throat is at 57.785 mm from the nozzle inlet. Fig. 1b shows the computational domain for the present problem. The present numerical analysis is performed at 2D symmetric condition. The nozzle geometry with grids are shown in Fig. 1c. The upper and downstream boundaries are at 40ht and 100ht from the throat of the nozzle. This domain size is found independent to the flow reflection from the boundaries. The nozzle pressure ratio (NPR) p01/pb, where p01 is inlet total pressure and pb is the back pressure, is oscillated sinusoidally as shown in Fig. 2. The back pressure is kept constant at 101.3 kPa. The oscillation of NPR consists of two processes namely-increasing process (NPR = 2.0–3.0) and decreasing process (NPR = 3.0–2.0). Four cases are considered by varying the frequency of oscillation. These are: Case A-50 Hz, Case B-100 Hz, Case C-150 Hz, and Case D-200 Hz. The corresponding time periods are 20 ms, 10 ms, 6.67 ms, and 5 ms (Fig. 2). The inlet total temperature is kept at 298.15 K. The initial and boundary conditions are presented as follows: at the inlet boundary, stagnation pressure and temperature were imposed as physical boundary conditions and the other variables were resulted from the numerical boundary treatment through Riemann invariant. Further, the values of k and x are specified as 1.22 m2/s2 and 61979.3 1/s, respectively at the inlet of the domain. The exit boundary was constrained with pressure outlet boundary

A.B.M. Toufique Hasan / International Journal of Heat and Fluid Flow 46 (2014) 70–83

Point A B C D E F G H I

73

Coordinates (mm) x y 0 0 0 -15.4 0 34.7 22.9 29.1 49.7 15.3 59.9 13.8 60.8 14.0 56.9 29.2 113.8 24.3

(a) Pressure Outlet Ambient Pressure Inlet

Pressure Outlet Symmetry

Pressure Inlet

(b)

(c)

Fig. 1. (a) Nozzle geometry (Hunter, 2004). (b) Computational domain with boundary conditions. (c) Close view of C–D nozzle with grids.

: Case-A

: Case-C

: Case-B

: Case-D

However, the differences in shock locations are within 0.6% in the computations with grid sizes of 76,000 (750  60) and 85,400 (900  60). Considering the computational time, a grid size of 76,000 (750  60) was found optimum and used for all the computations in the present study. In this case, the minimum normal grid spacing was reduced to give the value of wall y + 0.5. Further, the unsteady computations of nozzle flows under oscillating NPR are performed for various time step sizes, Dt. The locations of incident shock waves (IS) at various NPRs with different time step sizes are shown in Fig. 3b. It is found that the shock locations are dependent on the time step size. However, the differences in shock locations are within 0.2% in the computations with Dt = 106 s and Dt = 108 s. Thus, a time step size of 106 s was considered sufficient for this type of unsteady computation. Moreover, 300 internal iterations are performed at each time step to ensure the converge of the solution variables.

2.5

2.0

0

5

g sin rea dec

inc rea sin g

NPR

3.0

10

15

20

time period, T (ms) Fig. 2. Variation of NPR with time.

3. Computational validation condition. Non-slip and adiabatic wall conditions were applied at the solid boundary. The pressure at the wall was obtained from zero normal pressure gradients on the body surface. Moreover, the treatment of turbulence quantities at the wall is performed by the following equation (Wilcox, 2010):

k ¼ 0 and x ¼ 10

6m bðDy1 Þ2

ð6Þ

where m = ll/q, b is a turbulence model constant of value of 0.075 and Dy1 = distance of the first grid point away from the wall (0.90 lm in the present computation). Computational domain is discretized by structured mesh (Fig. 1b). Numerical computations are performed with different grid sizes to ensure the grid independent solution. The computed wall pressure distributions are shown in Fig. 3a. In this figure, the number of grids shown in the bracket is for diverging portion of the C–D nozzle and the total number of grids is shown outside the bracket. It was found that the magnitude of flow expansion and the location of shock wave are dependent on the grid sizes.

The performance of the present computational methods is verified against experimental data. A rigorous comparison between present computation and Hunter’s (2004) experiment in pressure ratios of 2.4 and 3.0 are implemented. Fig. 4 shows the experimental and numerical schlieren images (density gradient, rq) obtained from the present computation. The abscissa is normalized by the throat position (xt) measured away from the nozzle inlet. For both the cases of NPR = 2.4 (Fig. 4a) and 3.0 (Fig. 4b), the computed flow structures are captured reasonably with the experimental data in overexpanded nozzle flow conditions. The flow field comprises a well defined lambda shock structure and a free shock separation (FSS) that begins at the leading lambda shock and extends downstream. The locations of shocks and Mach stems are manually extracted from the numerical schlieren images in the computation. The computed leading lambda shocks located at x/xt  1.42 and 1.58 for NPR = 2.4 and 3.0, respectively. Moreover, the computationally obtained Mach stems located at x/xt  1.70 and 1.95 for NPR = 2.4 and 3.0, respectively. These locations are within maximum relative error of 3% of experimental data. The detail shock

A.B.M. Toufique Hasan / International Journal of Heat and Fluid Flow 46 (2014) 70–83

1.7

1 : 43,200 (400x50) : 55,600 (500x60) : 76,400 (750x60) : 85,400 (900x60)

1.6 1.5

p/p0

xIS /xt

0.6 0.4

1.4

g

0.8

: time step =10-4 : time step =10-6 : time step =10-8

i nc rea si n

74

1.3 0.2 1.2 0 0.6

0.8

1

1.2

1.4

1.6

1.8

2

1.8

2

2.2

2.4

2.6

x/xt

NPR

(a)

(b)

2.8

3

3.2

Fig. 3. (a) Distribution of wall pressure for different grid sizes at NPR = 3.0 and (b) effect of time step size on unsteady computations.

Experimental

(a) NPR = 2.4

Experimental

Present computation

Present computation

(b) NPR = 3.0 Fig. 4. Comparisons of experimental (Hunter, 2004) and present computational Schlieren images for flow through a converging–diverging nozzle.

Reflected Shock 50°

Mach Stem 71°

18.6 mm

Triple Point Incident Shock

(a)

(b)

Fig. 5. Shock schematic at NPR = 2.4 from. (a) Experiment (Hunter, 2004). (b) Present computation.

75

A.B.M. Toufique Hasan / International Journal of Heat and Fluid Flow 46 (2014) 70–83

0.25

wall shear stress (Pa)

1

p/p0

0.8 0.6 0.4 0.2

0.2 0.15 NPR 2.4 NPR 3.0

0.1 0.05 0

0

0.5

1

1.5

2

0

x/xt :Experiment-NPR 2.4 :Experiment-NPR 3.0

0.5

1

1.5

2

x/xt

:Computation-NPR 2.4 :Computation-NPR 3.0

(a)

(b)

Fig. 6. Distribution of (a) wall pressure and (b) wall shear stress for converging–diverging nozzle.

(f) NPR = 3.0

(e) NPR = 2.8

(g) NPR = 2.8

(d) NPR = 2.6

(h) NPR = 2.6

(c) NPR = 2.4

(i) NPR = 2.4

(b) NPR = 2.2

(j) NPR = 2.2

(a) NPR = 2.0

(k) NPR = 2.0

Fig. 7. Isocontours of density during increasing and decreasing processes of NPR oscillation for flow through a converging–diverging nozzle (Case-B, f = 100 Hz).

76

A.B.M. Toufique Hasan / International Journal of Heat and Fluid Flow 46 (2014) 70–83

structure for NPR = 2.4 is shown in Fig. 5 both from experiment and computation. Actually, the leading lambda shock comprises of an incident shock (IS) wave of oblique type which originates from the nozzle wall and a corresponding reflected shock (RS). These shock waves meet with Mach steam (MS) at the triple point as shown in the figure. Detail shock structures are well captured from computational results. Further, to quantify the properties of flow structure, the diameter of Mach stem and the two angles- one between the nozzle wall and incident shock and the other between the incidence and reflected shocks are considered as the comparison criteria. The computed diameter and angles are within relative error of 3.6% and 3.9%, respectively compared to experimental measurements. A comparison between the computed and experimental dimensionless wall static pressures plotted as a function of normalized distance is shown in Fig. 6a. Pressure is non-dimensionalized by chamber stagnation pressure (p01). The axial distribution of pressure is similar in both the cases of experiment and computation. However, the computational separation occurs at x/xt  1.42 which is only 0.7% different from the experimental results (x/xt = 1.41) for

NPR = 2.4. In case of NPR = 3.0, the variation is within 1.3%. There are minor differences between the computation and experimental results because of the assumption of 2D symmetric geometry even though the experiment was performed in the non-axisymmetric three-dimensional nozzle. The pressure peak found in Fig. 6a justifies the formed separation, its accompanied oblique shock, and the existence of adverse pressure gradient. Further, the distribution of wall shear stress is shown in Fig. 6b to quantify the location of separation points. The locations of negative wall shear stress confirmed the same locations of FSS for NPR = 2.4 and 3.0. Finally, it can be concluded that a satisfactory agreement exists between the present computation and the Hunter’s experiment.

4. Results and discussion Fig. 7 shows the isocontours of density during the successive increasing and decreasing modes of change of pressure ratio for oscillation frequency, f = 100 Hz (Case-B). The left side sequence presents the increasing process of NPR and the right side for the

(f) NPR = 3.0

(e) NPR = 2.8

(g) NPR = 2.8

(d) NPR = 2.6

(h) NPR = 2.6

(c) NPR = 2.4

(i) NPR = 2.4

(b) NPR = 2.2

(j) NPR = 2.2

(a) NPR = 2.0

(k) NPR = 2.0

Fig. 8. Isocontours of density during increasing and decreasing processes of NPR oscillation for flow through a converging–diverging nozzle (Case-D, f = 200 Hz).

77

A.B.M. Toufique Hasan / International Journal of Heat and Fluid Flow 46 (2014) 70–83

decreasing process of pressure ratio. At the beginning of increasing process of NPR = 2.0 (Fig. 7a), the incident shock wave and Mach stem are located at x/xt  1.28 and 1.51, respectively. The incident shock wave interacts with the boundary layer at nozzle upper wall and the corresponding free shock separation (FSS) can be seen in the figure. Moreover, this viscous separation interact with the

inviscid subsonic core and generates flow disturbances downstream of the Mach stem. With successive increasing of NPR, the flow structure moves downstream continuously as seen in Fig. 7b–e for NPR = 2.2, 2.4, 2.6, and 2.8, respectively. However, at higher NPRs of 2.6 and 2.8, the shock cells are well developed and visible downstream the Mach stem. At the end of the

(a) NPR = 2.1

(b) NPR = 2.2

(c) NPR = 2.4

(d) NPR = 2.6

(e) NPR = 2.8

(f) NPR = 2.9

Fig. 9. Sonic line inside the nozzle during increasing and decreasing processes of NPR oscillation (Case-B, f = 100 Hz).

1.7

1.6

i nc rea sin

g

1.6

2.2

2.4

2.6

2.8

3

in de

2.2 2.4 2.6 2.8

g

(b) f = 100Hz

1.4

1.3

1.3

1.2

1.2 2.2 2.4 2.6 2.8

3

3.2

1.8

3.2

g

inc rea si n

i

xIS /xt

de

cr ea sin

g

1.5

g s in ea r nc

3

g

(a) f = 50Hz

1.6

2

2

NPR

1.6

1.8

cr

ea s 1.8

3.2

NPR

de cre asi n

2

1.7

1.4

i

in as re c n

1.2

1.7

1.5

xIS /xt

1.4 1.3

1.3 1.8

g

1.5

xIS /xt

g 1.4

de cre asi n

xIS /xt

1.5

2

2.2 2.4 2.6 2.8

NPR

NPR

(c) f = 150Hz

(d) f = 200Hz

3

3.2

Fig. 10. Locations of incident shock waves during increasing and decreasing processes of NPR oscillation.

A.B.M. Toufique Hasan / International Journal of Heat and Fluid Flow 46 (2014) 70–83

increasing process of NPR = 3.0 (Fig. 7f), the incident shock and the Mach stem are located at x/xt  1.57 and 1.97, respectively. After that the NPR starts to decrease. However, the response of pressure decrease on the flow structure is not spontaneous. At NPR = 2.8 (Fig. 7g), the incident shock and the Mach stem are located at x/xt  1.56 and 1.96, respectively which are quite downstream to the corresponding locations of x/xt = 1.53 and 1.79 in increasing process (Fig. 7e). Then the flow structure is shifted upstream continuously at successive decreasing of NPR as shown in Fig. 7h–j. At the end of the decreasing process of NPR = 2.0 (Fig. 7k), the flow structure becomes almost the same as that in increasing process (Fig. 7a). Thus completes the cycle of pressure ratio oscillation. As seen from these figure sequences, there are differences between the locations of the incident shock waves (IS) and Mach stem (MS) even the nozzle at the same NPR. Results of isocontours of density for the case of f = 200 Hz (CaseD) are shown in Fig. 8. In this case, at the beginning of increasing process at NPR = 2.0 (Fig. 8a), the incident shock wave and Mach stem are located at x/xt  1.26 and 1.49, respectively. These locations are slightly upstream compared to the case of f = 100 Hz (Fig. 7a). Moreover, at next NPR of 2.2 (Fig. 8b) in the successive increasing process, these locations are moved further upstream to x/xt  1.22 and 1.44 though the flow is ideally at more expanded conditions compared to NPR = 2.0 (Fig. 8a). Thus it can be seen that

at higher oscillating condition, the response of initial pressure ratio increment on flow structure is seen to have a lag compared to lower frequencies (for example in Fig. 7). After that, the flow responses spontaneously with the increment of pressure ratio as seen in Fig. 8c–e for NPR = 2.4, 2.6 and 2.8, respectively. At NPR = 3.0 (Fig. 8f), the shock and Mach stem are located at x/xt  1.58 and 1.99, respectively. Thus ends the increasing process. Then, the NPR starts to decrease and at the initial NPR = 2.8 (Fig. 8g), the locations of shock and Mach stem are at x/xt  1.59 and 1.98. The corresponding locations are at x/xt  1.54 and 1.90 in increasing process (Fig. 8e). With further decrease in pressure ratio, the flow structure is shifted upstream as seen in Fig. 8h–j. Finally, the flow resemble the same flow structure at NPR = 2.0 (Fig. 8k) as compared to Fig. 8a. Thus, the irreversibility in the location of flow structure is confirmed from both the flow cases. The flow structure for the cases of f = 50 Hz and 100 Hz also correspond to the similar irreversible behavior during the oscillation of NPR. These sequential figures are not shown here for brevity. However, the specific locations of incident shock wave and Mach stem will be addressed in the following discussion. Now the distribution of sonic (M = 1) lines inside the nozzle under the imposed NPR oscillation will be discussed. Results are shown for oscillation frequency of 100 Hz as shown in Fig. 9. The solid line is for increasing process and the dotted line is for

2

2

1.9

1.9 1.8 1.7

1.6

1.6

1.5

1.5

1.4

1.4 2.2 2.4 2.6 2.8

3

1.8

3.2

(a) f = 50Hz

(b) f = 100Hz

1.9

1.9

in cr ea s

1.7

1.6

1.6

1.5

1.5

in

1.8

3.2

g

g in

in de cre as

3

g

2

g

2.1

2

1.7

2.2 2.4 2.6 2.8

NPR

2.1

1.8

2

NPR

g

in cr ea s

2

in

sin ea cr

in

1.8

xMS /xt

g

cr ea si n

i

in as re nc

de

1.7

g

xMS /xt

in as re c de

1.8

g

2.1

xMS /xt

xMS /xt

2.1

de cre as

78

1.4

1.4 1.8

2

2.2 2.4 2.6 2.8

3

3.2

1.8

2

2.2 2.4 2.6 2.8

NPR

NPR

(c) f = 150Hz

(d) f = 200Hz

3

3.2

Fig. 11. Locations of Mach stem xMS during increasing and decreasing processes of NPR oscillation.

Table 1 Percentage deviation of the location of Mach stem, DxMS/xt in decreasing process during NPR oscillation. NPR Oscillation frequency

50 Hz 100 Hz 150 Hz 200 Hz

2.0

2.1

2.2

2.3

2.4

2.5

2.6

2.7

2.8

2.9

3.0

– – – –

1.5 1.6 2.8 6.2

3.5 4.2 5.7 12.5

4.3 5.5 7.3 13.9

4.6 5.2 7.5 13.4

4.1 5.0 7.4 11.4

3.5 4.8 7.3 9.7

3.1 4.2 6.7 7.7

2.1 3.6 5.3 5.3

1.2 2.0 3.0 2.5

– – – –

79

A.B.M. Toufique Hasan / International Journal of Heat and Fluid Flow 46 (2014) 70–83

decreasing process. For all the NPRs, sonic lines are observed at the throat which are concave in shape while looking from left side. Other sonic line is seen to generate from the point of flow separation from the wall. And another at the location of Mach stem which is convex in shape while looking from the same direction. These two sonic lines meet at the triple point and then distribute downstream as slip sonic line. There is no significant difference in the sonic line structure is observed during the NPR oscillation. However, difference in the location of sonic lines are observed (except at the throat position) at the same NPRs during the oscillation. At NPR = 2.1 (Fig. 9a), the sonic lines are located at x/xt  1.48 and 1.52 for increasing and decreasing processes, respectively. At higher NPRs, the differences are successively increases (as shown in Fig. 9b and c) and it is maximum at NPR = 2.6 (Fig. 9d). At this NPR, the locations of sonic lines are at x/xt  1.75 and 1.90 for increasing and decreasing processes, respectively. After that the differences are decreases with increase of NPRs as shown in Fig. 9e and f for NPR = 2.8 and 2.9, respectively. Other lines – for example sonic lines from separation point and triple point, the differences are observed in the vertical direction. However, the differences in locations during the oscillation are in the same manner as discussed earlier. Fig. 10 shows the locations of incident shock waves (IS) which are developed from the nozzle wall at overexpanded flow conditions during the oscillation of NPR. Thick line with open circle presents the locations of successive increasing process and the thin line with closed circle presents the corresponding locations at successive decreasing process. For all the cases of oscillation frequencies, the flow inside the nozzle experiences a lag effect at initial NPRs during the increasing process. The upstream movement of incident shock waves at initial NPRs from the corresponding locations at NPR = 2.0 is considered to be due to this effect. This effect can be seen up to NPR = 2.10, 2.20, 2.20 and 2.30 for the cases of

f = 50 Hz (Fig. 10a), 100 Hz (Fig. 10b), 150 Hz (Fig. 10c), and 200 Hz (Fig. 10d), respectively. After this initial lag, the shock waves are generated downstream spontaneously in response to the increasing process of NPR. After the end of increasing process at NPR = 3.0, the flow experiences a lead effect during the initial NPRs of decreasing process. This lead effect causes to generate the shocks further downstream though the flows are ideally to be in less expanded conditions corresponding to NPR = 3.0. The lead effect exists up to NPR = 2.90, 2.80, 2.80 and 2.70 for the cases of f = 50 Hz (Fig. 10a), 100 Hz (Fig. 10b), 150 Hz (Fig. 10c), and 200 Hz (Fig. 10d), respectively. Then, the response of shock system to the decrement of NPRs becomes spontaneous. And finally, the flow structure merges with the state of NPR = 2.0 from where the NPR oscillation of next cycle starts again. These lag and lead effects cause the different shock locations at the same NPRs during the oscillation of NPR. The shocks are at downstream locations in decreasing process compared to the increment during pressure oscillation. Adding all the successive locations, a hysteresis loop can be found for all the flow cases as shown in Fig. 10. The size of the hysteresis loop is seen to increase at higher oscillation frequencies. The computed maximum downstream movement of shock locations are 1.6% at NPR = 2.4, 6.9% at NPR = 2.3, 8.0% at NPR = 2.3, and 14% at NPR = 2.3 for f = 50 Hz, 100 Hz, 150 Hz, and 200 Hz, respectively. The variations of the locations of Mach stem (MS) with NPR oscillation are shown in Fig. 11. In this case, the lag and lead effects are seen to be slightly smaller as compared to Fig. 10. However, a hysteretic behavior in the locations of Mach stem are well captured for all the flow cases of f = 50 Hz (Fig. 11a), 100 Hz (Fig. 11b), 150 Hz (Fig. 11c), and 200 Hz (Fig. 11d). The Mach stems are located at downstream locations in decreasing process compared to the increasing process at the same NPR during the oscillation. The percentage deviation of the location of Mach stem, DxMS/xt

2.2

2.2

2

2 1.8

si n rea c de

1.6 1.4

g

in

as cre

ing

pMS /pb

pMS /pb

1.8

1.6 1.4

1.2

1.2

1

1

0.8 1.8

2

2.2

2.4

2.6

2.8

3

0.8 1.8

3.2

d

s in rea ec

2

g

2.2

NPR

2.2

2

2

1.4

1.8

ng i

i as re c n

ng

pMS /pb

pMS /pb

1.6

1.6 1.4

1.2

1.2

1

1

0.8 1.8

2

2.2

2.4

2.6

2.6

2.8

3

3.2

(b) f = 100Hz

2.2

ii as re c de

2.4

g

NPR

(a) f = 50Hz

1.8

sin rea c in

2.8

3

3.2

0.8 1.8

in as re c de

g

in

2

2.2

2.4

2.6

2.8

NPR

NPR

(c) f = 150Hz

(d) f = 200Hz

sin ea cr

3

g

3.2

Fig. 12. Strength of Mach stem pMS during increasing and decreasing processes of NPR oscillation.

80

A.B.M. Toufique Hasan / International Journal of Heat and Fluid Flow 46 (2014) 70–83

during the NPR oscillation is shown in Table 1. DxMS/xt is calculated as follows:

% DxMS =xt ¼

DpMS/pb during the NPR oscillation. The definition of DpMS/pb is kept similar as defined in the preceding section and shown as below

ðxMS =xt Þdecreasing process  ðxMS =xt Þincreasing process  100% ðxMS =xt Þincreasing process

% DpMS =pb ¼

ð7Þ

ðpMS =pb Þdecreasing process  ðpMS =pb Þincreasing process  100% ðpMS =pb Þincreasing process ð8Þ

Results show that the percentage deviation of Mach stem location is increased from zero, reaches a peak and finally decreases to zero during the NPR oscillation. The peaks of deviation of Mach stem locations are 4.6%, 5.5%, 7.5%, and 13.9% for the cases of f = 50 Hz, 100 Hz, 150 Hz, and 200 Hz, respectively. The peak deviations are occurred at NPR = 2.4 for f = 50 Hz and at NPR = 2.3 for all the remaining cases. Fig. 12 shows the effect of NPR oscillation on the strength of Mach stem. The strength is defined as the ratio of the pressure at the Mach stem (pMS) and that at the back flow condition (pb). A hysteretic behavior in the strength of Mach stem are also observed for all the flow cases of f = 50 Hz (Fig. 12a), 100 Hz (Fig. 12b), 150 Hz (Fig. 12c), and 200 Hz (Fig. 12d). Mach stem becomes strengthen in the decreasing process. It is found that the Mach stem strength showed the same hysteresis behavior with higher magnitude compared to other parameters as discussed earlier. To better understand the trend of hysteretic characteristics, Table 2 presents the percentage deviation of the strength of Mach stem,

In this consideration, the deviation shows the similar pattern as in the locations of Mach stem. However, the magnitudes of peak deviations are 4.4%, 10.0%, 19.2%, and 30.5% for the cases of f = 50 Hz, 100 Hz, 150 Hz, and 200 Hz, respectively. The distribution of the turbulent viscosity, lt inside the nozzle flow field is shown in Fig. 13 to describe the role of shear stress transport (SST) of presently implemented turbulence model. Results are shown for only case of 200 Hz oscillation for brevity. It is observed that the turbulent viscosity is induced only at the zone of shock-boundary layer interaction (SBLI). Thus the turbulent viscosity evolves just after the location of free-shock separation and after that the magnitude of the viscosity is gradually increased. It is observed that the ratios of turbulent viscosity to laminar viscosity, lt/ll are pretty high in the region of SBLI. The values of viscosity ratios are found to be up to 3000 in the region of interest. Since the distribution of turbulent viscosity follows characteristics of shock structure, the distribution of viscosity is quite different even

Table 2 Percentage deviation of the strength of Mach stem, DpMS/pb in decreasing process during NPR oscillation.

50 Hz 100 Hz 150 Hz 200 Hz

2.0

2.1

2.2

2.3

2.4

2.5

2.6

2.7

2.8

2.9

3.0

– – – –

2.4 4.0 8.3 11.7

2.8 7.8 14.3 19.1

4.2 10.0 19.2 30.5

4.4 9.0 18.0 28.8

3.7 8.1 15.8 27.9

3.5 7.6 14.8 20.1

2.2 7.0 9.3 16.4

1.8 6.0 8.4 13.0

1.5 3.5 5.0 3.7

– – – –

y/xt

0.4

0.02

0.01

0.04

0.03

5

Oscillation frequency

0.0

NPR

0.2 0 1

1.5

2

2.5

x/xt

(c) NPR = 3.0 0.4 0.02 0.03

0.4 0.04

0.01

0.05

y/xt

y/xt

0.01

0.2 0

0.03

0.04

0.05

0 1

1.5

2

1

2.5

1.5

2

2.5

x/xt

x/xt

(d) NPR = 2.6

(b) NPR = 2.6 0.4

0.4 0.01 0.02 0.03

0.04

0.01

0.05

y/xt

y/xt

0.02

0.2

0.2

0.02 0.03 0.04

0.05

0.2 0

0 1

1.5

2

2.5

1

1.5

(a) NPR = 2.4 Level 1 TurbulentViscosity: 0.00 (kg/m-s) t

2

2.5

x/xt

x/xt

(e) NPR = 2.4 2 0.01

3 0.02

4 0.03

5 0.04

6 0.05

Fig. 13. Distribution of turbulent viscosity lt inside the nozzle during increasing and decreasing processes of NPR oscillation (Case-D, f = 200 Hz).

81

A.B.M. Toufique Hasan / International Journal of Heat and Fluid Flow 46 (2014) 70–83

(a) 0.8

0.8

sin ea cr e d

0.7

g in

0.6

0.5 1.8

2

2.2

2.4

2.6

si n ea cr e d

F / Fi

F / Fi

0.7

ing as cre

g sin rea i nc

0.6

2.8

3

3.2

0.5 1.8

2

2.2

2.4

2.6

NPR

NPR

(b)

(c)

2.8

g

3

3.2

Fig. 14. (a) Procedure for calculating thrust (Papamoschou et al., 2009), and thrust coefficient during increasing and decreasing processes of NPR oscillation for the cases of (b) f = 50 Hz, and (c) f = 200 Hz.

the nozzle at the same NPR during the decreasing and increasing processes of oscillation. The thrust of a nozzle equals the integral of the axial components of the forces acting on the internal and external walls of the nozzle and supply reservoir. The procedure for calculating the thrust is shown in Fig. 14a (Papamoschou et al., 2009). The thrust coefficient for the actual nozzle is given by the following equation (Papamoschou et al., 2009):

F 1 ¼1þ _ e Fi mU

Z

e

ðp  pa Þ sin hds

ð9Þ

d

where F is the actual thrust and Fi is the ideal thrust for perfectly expanded flow condition. p is the internal pressure and pa is the external pressure at static condition. h is the local angle of the nozzle surface. The mass flow rate is based on the sonic conditions in the actual nozzle,

_ ¼ At c m



p01 2 a01 c þ 1

ccþ1 1 ð10Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffi where a01 ¼ cRT 01 is the speed of sound at the reservoir. The fully-expanded velocity is given by the following relation

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 U e ¼ a0 ð1  NPRð1cÞ=c Þ c1

ð11Þ

Integration of pressure started at the point ‘‘d’’ where p = pa and ended at the nozzle exit. The variation of thrust coefficient during the increasing and decreasing process of the NPR oscillation is shown in Fig. 14b and c for the cases of 50 Hz and 200 Hz, respectively. It is found that the thrust coefficient increases with increase of NPR for both the cases. However, there are no significant differences are observed in the thrust performance during the increasing and decreasing processes of NPR oscillation. For the case of 50 Hz oscillation (Fig. 14b), the maximum difference is within the range of 0.3%. In case of 200 Hz oscillation, the largest deviation is found in the order of 1.3% at NPR = 2.5 (Fig. 14c). Other two cases show the similar nature as in the case of 200 Hz. Thus, it can be

concluded that the overall nozzle thrust performance does not show the same order of irreversibility as in the cases of locations of shock structures. 5. Conclusions A numerical simulation is performed to investigate the effect of the oscillation of nozzle pressure ratio (NPR) on the shock structure inside a two-dimensional, non-axisymmetric converging–diverging (C–D) supersonic nozzle. Overexpanded flow conditions are considered which is dominated by the free shock separation (FSS). The effect of oscillation frequency is also explored in this study. Computational results are well validated with the available experimental data. The results of the present study can be summarized as follows: (i) The flow structure exhibits an irreversible characteristics during the NPR oscillation. The locations of free shock separation, Mach stem, triple points and the distribution of sonic lines are significantly different in increasing and decreasing processes even at the same NPR during the oscillation of nozzle pressure ratio. (ii) The successive locations of the shock structure during increasing and decreasing processes produces a hysteresis loop in a complete cycle of NPR oscillation. (iii) The maximum differences in the locations of free shock separation are 1.6%, 6.9%, 8.0% and 14% in between increasing and decreasing processes for f = 50 Hz, 100 Hz, 150 Hz, and 200 Hz, respectively. (iv) The corresponding differences in the locations of Mach stem are 4.6%, 5.5%, 7.5% and 14% for f = 50 Hz, 100 Hz, 150 Hz, and 200 Hz, respectively. (v) The strength of Mach stem also presents the same hysteretic characteristics during the NPR oscillation. And the strengths are 4.4%, 10%, 19.2%, and 30.5% higher in decreasing process compared to increasing process at the same NPR. (vi) The computed results shows that the hysteretic characteristics are increased at higher oscillation frequencies.

82

A.B.M. Toufique Hasan / International Journal of Heat and Fluid Flow 46 (2014) 70–83

(vii) However, the nozzle thrust performance is not affected significantly due to the NPR oscillation in comparison to the locations of flow structures.

where Re0 is the Reynolds number based on reference inlet flow conditions, Gk is the generation of turbulent kinetic energy and J is the Jacobian of transformation, and it is given by

J ¼ jxn yg  yn xg j Acknowledgements The present work has been carried out with computational resource support from Higher Education Quality Enhancement Project (HEQEP), AIF (2nd Round)-Sub-Project CP 2099, UGC, MoE, Government of Bangladesh (Contract No. 28/2012).

Other details of the turbulence model can be found in Menter (1994), Wilcox (2010). In numerical simulations, the contravariant velocities U and V are used as the velocity components in the generalized coordinates (n, g), which are related to the original velocities in Cartesian coordinates (u and v):

U ¼ ðu  xs Þnx þ ðv  ys Þny ¼ nt þ unx þ v ny

Appendix A

V ¼ ðu  xs Þgx þ ðv  ys Þgy ¼ gt þ ugx þ v gy

The governing equations of the present problem are first normalized with the reference values at the inlet flow conditions and then mapped into computational space of curvilinear coordinate system (n, g) using the generalized transformation:

n ¼ nðx; y; tÞ g ¼ gðx; y; tÞ

s¼t

ðA:1Þ

Applying the transformation the governing equations can be written in the flux vector form:

b @E b @ Fb b @ Fb @U 1 @R þ þ ¼ þ @ s @n @ g Re0 @n @ g

! b þH

ðA:2Þ

q

3

6 6 qu 6 6 6 qv b U ¼ J6 6 6 E 6 6 qk 4

7 7 7 7 7 7 7 7 7 7 5

2

3

2

nt ¼ xs nx  ys ny gt ¼ xs gx  ys gy

qU 7 6 6 quU þ nx p 7 7 6 7 6 6 qv U þ ny p 7 b 7 6 E ¼ J6 7 6 UðE þ pÞ 7 7 6 6 qkU 7 5 4

qx

qxU 0

qV 7 6 6 quV þ gx p 7 7 6 7 6 6 qv V þ g y p 7 b 7 6 F ¼ J6 7 6 VðE þ pÞ 7 7 6 6 qkV 7 5 4

The contravariant velocity components U and V are in directions normal to the constant n and g surfaces, respectively. The quantity (xs and ys) is the velocity of any points on the ‘‘grid’’ in an initial frame. In the present work, the body is not in motion and these velocities are zero. The viscous flux terms with the shear stresses in the transformed coordinates are:

0

6 6 6 6 6 6 b H ¼ J6 6 6 6 6 4

ðA:8Þ

nx ¼

ðA:9Þ

1 y ; J g

1 n y ¼  xg ; J

1 J

1 J

gx ¼  yn ; gy ¼ xn

ðA:10Þ

The auxiliary functions are:

leff @T ðc  1ÞPr @x   leff @a2 @a2 ¼ sxx u þ sxy v þ nx þ gx ðc  1ÞPr @n @g

a ¼ sxx u þ sxy v þ

leff @T ðc  1ÞPr @y   leff @a2 @a2 ¼ syx u þ syy v þ ny þ gy ðc  1ÞPr @n @g

ðA:11Þ

b ¼ syx u þ syy v þ

ðA:12Þ

Pr is the Prandtl number, a is the speed of sound, and a2 = cRT. c is the specific heat ratio and a value of 1.40 has been used in the present study.

3

0 0 0 0

Gk  qb kxRe0 qa lt Gk

sxy ¼ leff ðun ny þ ug gy þ v n nx þ v g gx Þ

The quantities nx, ny, gx, gy are called the metrics of transformation, defined as follows:

3

7 6 7 6 gx sxx þ gy sxy 7 6 7 6 7 6 gx sxy þ gy syy 7 6 7 6 b S ¼ J6 7 g a þ g b x y 7 6  7 6  6 l þ lt ðg k þ g k Þ 7 7 6 x y x y l r k 7 6  5 4 lt ll þ rx ðgx xx þ gy xy Þ 2

ðA:7Þ

2 3

7 6 7 6 nx sxx þ ny sxy 7 6 7 6 7 6 n s þ n s xy yy x y 7 6 7 b ¼ J6 R 7 6 n a þ n b x y 7 6   7 6 6 l þ lt ðn k þ n k Þ 7 7 6 x x y y l rk 7 6  5 4 lt ll þ rx ðnx xx þ ny xy Þ 2

sxx ¼ leff ½2ðun nx þ ug gx Þ  ðv n ny þ v g gy Þ

syy ¼ leff ½2ðv n ny þ v g gy Þ  ðun nx þ ug gx Þ

qxV

3

ðA:6Þ

2 3

3

2

ðA:5Þ

where

where

2

ðA:4Þ

1  qbx2 Re0  2ð1  F 1 Þqrx;2 xRe ðkx þ ky Þðxx þ xy Þ 0

7 7 7 7 7 7 7 7 7 7 7 5

ðA:3Þ

References Balabel, A., Hageb, A.M., Nasr, M., El-Behery, S.M., 2011. Assessment of turbulence modeling for gas flow in two-dimensional convergent–divergent rocket nozzle. Appl. Math. Model. 35, 3408–3422. Barth, T.J., Jespersen, D., 1989. The design and application of upwind schemes on unstructured meshes. In: Proceedings of the Twenty-Seventh AIAA Aerospace Science Meeting, Reno, Nevada, AIAA-89-0366. Ben-Dor, G., Vasiliev, E.I., Elperin, T., Chpoun, A., 2001. Hysteresis phenomena in the interaction process of conical shock waves: experimental and numerical investigations. J. Fluid Mech. 448, 147–174. Ben-Dor, G., Elperin, T., Vasiliev, E.I., 2003. Flow-Mach-number-induced hysteresis phenomena in the interaction of conical shock waves – a numerical investigation. J. Fluid Mech. 496, 335–354.

A.B.M. Toufique Hasan / International Journal of Heat and Fluid Flow 46 (2014) 70–83 Bruce, P.J.K., Babinsky, H., 2008. Unsteady shock wave dynamics. J. Fluid Mech. 603, 463–473. Bur, R., Benay, R., Galli, A., Berthouze, P., 2006. Experimental and numerical study of forced shock-wave oscillations in a transonic channel. Aerosp. Sci. Technol. 10, 265–278. Burstchell, Y., Zeitoun, D., Chpoun, A., Ben-Dor, G., 2001. Conical shock interactions in steady flows: Numerical study. In: Proceedings of the Twenty Third International Symposium on Shock Waves, Texas. Frey, M., Hagemann, G., 1999. Flow separation and side-loads in rocket nozzles. In: Proceeding of Thirty-Fifth AIAA/ASME/ASEE Joint Propulsion Conference and Exhibit, Los Angeles, CA AIAA paper-99-2815. Gatski, T.B., Bonnet, J.-B., 2013. Compressibility, Turbulence, and High Speed Flow. Elsevier. Hadjadj, A., Onofri, M., 2009. Nozzle flow separation. Shock Waves 19, 163–169. Hadjadj, A., Kudryavtsev, A.N., Ivanov, M.S., 2004. Numerical investigation of shockreflection phenomena in overexpanded supersonic jets. AIAA J. 42, 570–577. Hunter, C.A., 2004. Experimental investigation of separated nozzle flows. J. Propul. Power 20, 527–532. Johnson, A., Papamoschou, D., 2010. Instability of shock-induced nozzle flow separation. Phys. Fluids 22, 016102. Lijo, V., Kim, H.D., Setoguchi, T., Matsuo, S., 2010. Numerical simulation of transient flows in a rocket propulsion nozzle. Int. J. Heat Fluid Flow 31, 409–417. McIntosh, A., Rylands, S., 1996. A model of heat transfer in Rijke tube burners. Combust. Sci. Technol. 113, 273–289. Menter, F.R., 1994. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J. 32, 269–289.

83

Olson, B.J., Lele, S.K., 2013. A mechanism for unsteady separation in over-expanded nozzle flow. Phys. Fluids 25, 110809. Onofri, M., Nasuti, F., 2001. Theoretical considerations on shock reflections and their implications on the evaluation of air intake performance. Shock Waves 11, 151– 156. Papamoschou, D., Zill, A., Johnson, A., 2009. Supersonic flow separation in planar nozzles. Shock Waves 19, 171–183. Roe, P., 1981. Approximate Riemann solvers, parameters vector, and difference schemes. J. Comput. Phys. 43, 357–372. Roe, P., 1986. Characteristics based schemes for the Euler equation. Annu. Rev. Fluid Mech. 18, 337–365. Semlitsch, B., Mihaescu, M., Gutmark, E.J., Fuchs, L., 2013. Flow structure generation by multiple jets in supersonic cross-flow. In: Proceedings of 4th International Conference on Jets, Wakes and Separated Flows, ICJWSF2013, Nagoya, Japan. Shimshi, E., Ben-Dor, G., Levy, A., 2011. Viscous simulation of shock reflection hysteresis in ideal and tapered overexpanded planar nozzle. Shock Waves 21, 205–214. Strasse, W., 2011. Towards the optimization of a pulsatile three-stream coaxial airblast injector. Int. J. Multiph. Flow 37, 831–844. Wilcox, C., 1986. Multiscale model for turbulent flows. In: Proceedings of the Twenty-Fourth AIAA Aerospace Science Meeting, Reno, Nevada, USA. Wilcox, C., 2010. Turbulence Modeling for CFD. DCW Industries Inc. Xiao, Q., Tsai, H.M., Papamoschou, D., 2009. Experimental and numerical study of jet mixing from a shock containing nozzle. J. Propul. Power 25, 688–696.