Prediction of the onset of supersonic inlet buzz

JID:AESCTE AID:105523 /FLA

[m5G; v1.261; Prn:8/11/2019; 10:22] P.1 (1-10)

Aerospace Science and Technology ••• (••••) ••••••

Contents lists available at ScienceDirect

Aerospace Science and Technology www.elsevier.com/locate/aescte

Prediction of the onset of supersonic inlet buzz Jun Yamamoto a,1 , Yoimi Kojima a,1 , Masaharu Kameda a,∗ , Yasushi Watanabe b , Atsushi Hashimoto b , Takashi Aoyama b a b

Department of Mechanical Systems Engineering, Tokyo University of Agriculture and Technology, Koganei, Tokyo 184-8588, Japan Aeronautical Technology Directorate, Japan Aerospace Exploration Agency, Chofu, Tokyo 182-8522, Japan

a r t i c l e

i n f o

Article history: Received 1 July 2019 Received in revised form 25 October 2019 Accepted 26 October 2019 Available online xxxx Keywords: Supersonic flow Inlet Computational fluid dynamics Wind tunnel test Delay differential equation

a b s t r a c t A computational fluid dynamics (CFD) simulation was conducted about the onset of supersonic inlet buzz. Then, a novel mathematical model for the onset was proposed with the aid of the CFD results. The simulation was conducted using a common approach in industrial applications, in which unsteady Reynolds-averaged Navier-Stokes equations with millions of computational grids were used. The numerical results were validated with the experimental data obtained from a wind tunnel test, in which the flow around an external-compression inlet model at a Mach number of 2.0 was tested. The inlet model consists of a wedge for the generation of an oblique shock wave, a subsonic diffuser with a rectangular cross-section for deceleration of the captured flow, and a flow plug for constriction of the exit of the diffuser. Both the numerical and experimental results show that the onset of buzz correlated well with the position of the shock wave triple point (STP), which means that the onset is closely related to the mass flux through the diffuser. The critical position of the STP for the occurrence of buzz was well predicted by this simulation. Therefore, the balance of the mass in the rear part of the diffuser was modeled. The modeling indicates that this system can be represented using a delay differential equation containing a delayed negative feedback due to the pressure waves propagating in the diffuser. The presence of this feedback mechanism is the origin of the self-sustained oscillation of the buzz. The analysis shows that the onset and frequency of the inlet buzz are well predicted using the derived equation. In general, the onset of buzz is characterized using three parameters: the sensitivities of mass flow rates into/out of the part to its mass and the time lag of changes between the inflow rate and the mass. © 2019 Elsevier Masson SAS. All rights reserved.

1. Introduction

The return of supersonic transport (SST) for passengers is envisaged owing to recent vigorous research projects, such as the Drop test for the Simplified Evaluation of Non-symmetrically Distributed sonic boom (D-SEND) in Japan [1] and the Quiet Supersonic Transport (QueSST) low-boom flight demonstrator in the USA [2]. These projects assume a cruise Mach number of approximately 2.0. In order to re-launch supersonic passenger flights, it is necessary to overcome two serious problems, an unacceptable level of noise due to sonic booms and insufficient cruise range due to poor fuel consumption, which remained when the Concorde was developed in the 1960s.

* 1

Corresponding author. E-mail address: [email protected] (M. Kameda). J. Yamamoto and Y. Kojima contributed equally to this work.

https://doi.org/10.1016/j.ast.2019.105523 1270-9638/© 2019 Elsevier Masson SAS. All rights reserved.

The supersonic inlet is an important device to arrange the air supplied to an engine of an SST. The inlet for the aircraft with a cruise Mach number of 2.0 is generally external compression system, which consists of a wedge (or a cone) followed by a subsonic diffuser, as shown in Fig. 1. The wedge is used to generate an envelope of shock waves (or Mach waves) by which the air is efficiently decelerated to subsonic speed with low energy loss. The methodology of inlet design has been developed over the past 70 years. The basis for this design was established from 1950s to 1960s [3], followed by recent activity to update the computational design tool [4], the new inlet concepts for improved supersonic performance [5], the proposal of advanced inlet design to reduce the noise associated with sonic booms [6], and an inlet designed for a rocket based combined cycle engine [7]. The inlet must be designed to ensure that the flow in the diffuser is stable over a wide range of engine operation. The stability of the flow, however, deteriorates in subcritical operation, in which the mass flow through the diffuser is below the value at its optimal operating condition. This remarkable instability of the flow

JID:AESCTE

2

AID:105523 /FLA

[m5G; v1.261; Prn:8/11/2019; 10:22] P.2 (1-10)

J. Yamamoto et al. / Aerospace Science and Technology ••• (••••) ••••••

Fig. 1. Schematic diagram of supersonic inlet buzz.

with self-sustained oscillation of shock waves is known as supersonic inlet buzz. Inlet buzz has been investigated since the 1950s, not only from the perspective of inlet design but also as an intrinsic interesting phenomenon in fluid dynamics [3]. In the practical inlet deign, the position of cowl lip (see Fig. 1) is determined the oblique shock from the wedge to fall on the lip at the maximum Mach number of the aircraft plus 0.2 to avoid the buzz [8]. This design concept inevitably increases the inlet spillage drag. Fuel consumption due to the spillage drag can be suppressed by the reduction of this margin for Mach number if the onset of buzz is accurately predicted. Rapid progress in computational fluid dynamics (CFD) stimulated us to attempt quantitative prediction of such a complicated unsteady flow by CFD. Pioneering work in the numerical simulation of inlet buzz flow was conducted by Newsome [9] for an axisymmetric external-compression inlet in the freestream with a Mach number of 2.0. He used the historical MacCormack scheme to solve the Reynolds-averaged Navier-Stokes (RANS) equations. Although he did not complete quantitative analysis of the flow, he demonstrated that the self-sustained oscillation of shock waves could be simulated by CFD. Lu and Jain [10] then simulated another axisymmetric inlet configuration at a freestream Mach number of 2.0. They used the Osher-Chakravarthy upwind total variation diminishing (TVD) scheme to solve the RANS equations with the Baldwin-Lomax (BL) turbulence model. They showed that the time history of pressure fluctuation and its dominant frequency in the inlet were fairly close to the experimental value obtained by Dailey [11]. More recently, Trapier et al. [12] simulated the buzz flow of a rectangular, mixed-compression inlet at a freestream Mach number of 1.8. They used the delayed detached-eddy simulation (DDES) method, which is a hybrid RANS/LES (large-eddy simulation) approach, and the Spalart-Allmaras (SA) model as the turbulence model in the RANS simulation and showed that the frequency spectrum of pressure oscillation obtained by their DDES simulation agreed quantitatively well with their wind tunnel test results [13]. They also conducted an unsteady RANS (URANS) simulation and found that the primary dominant frequency of the pressure oscillation could be captured even if the URANS approach was used. From the viewpoint of inlet design, it is crucial to accurately predict the onset of inlet buzz. However, the predictability of the onset of inlet buzz by CFD has not been sufficiently studied, while the self-excited oscillation of shock waves associated with unsteady flow in the inlet [14–23] and the detection technique for the onset [24,25] have been well analyzed. Trapier et al. [12] suggested that it is difficult to determine the criterion for buzz onset through their numerical simulation. We consider that the difficulty in the prediction of buzz occurrence does not come mainly from the selection of governing equations or insufficient accuracy in CFD simulation. Rather, it is probably due to its high sensitivity to many factors, such as the complicated configuration of the inlet (e.g., [13,21,26,27]). The exact mechanism of self-sustained oscillation associated with inlet buzz remains unclear. Historically, two types of mechanisms, the Ferri criterion [28] and the Dailey criterion [11], were

reported to induce buzz [3]. The Ferri criterion is the occurrence of small-amplitude shock oscillation (little buzz) when the vortex sheet from the intersection point of normal and oblique shocks moves across the cowl lip. The Dailey criterion is the occurrence of large-amplitude shock oscillation (big buzz) with the unsteady flow triggered by the flow separation in the duct. Recently, however, based on their wind tunnel experiment, Chen et al. [21] claimed that the little buzz and big buzz may share a common origin, which is different from the prevailing understanding. Therefore, in order to elucidate the predominant factor characterizing the onset of inlet buzz, it is useful to examine wind tunnel experiments and numerical simulations for a simplified inlet model. In the present study, we propose a novel mathematical model for the onset of supersonic inlet buzz with the aid of the computational fluid dynamics simulation. In order to accurately predict the occurrence of inlet buzz within the range of computational resources in industrial applications, we used URANS-based numerical calculation. We simulated the flow around the inlet equivalent to the model in the wind tunnel experiment conducted by the Japan Aerospace Exploration Agency (JAXA). In order to prevent difficulty due to the complexity of the configuration of the inlet, we used a simpler model than the model used in common research into inlet buzz flow. We evaluated the feasibility of quantitative prediction through comparison with this experiment result. Finally, we show that the onset of self-sustained oscillation is modeled well using a delay differential equation. 2. Wind tunnel test We analyze the feasibility of our numerical approach using a wind tunnel test for characterization of supersonic inlet buzz conducted at JAXA. A brief description of this test follows. 2.1. Model Fig. 2 shows (a) a schematic diagram of the inlet model and (b) a photograph of the model installed in the test section of the wind tunnel. The main dimensions of the model are listed in Table 1. This model is an external-compression inlet consisting of a wedge for generation of an oblique shock wave, a subsonic diffuser with a rectangular cross-section for deceleration of the captured flow, and a flow plug for constriction of the exit of the diffuser. The flow is designed that the uniform flow with M = 2.0 is firstly decelerated by the oblique shock wave to M = 1.64, then further decelerated to M = 0.66 by the terminal shock wave in front of the entrance of diffuser, and finally M = 0.43 at the exit. The wedge is installed separately from the diffuser in order to prevent the boundary layer on the wedge surface from flowing into the diffuser, which is a major source of complexity in supersonic inlet buzz [22]. The mounting angle of the upper wall of the diffuser is variable in order to alter the performance of the diffuser. Transparent optical windows are embedded in the upstream part of the side wall of the diffuser. Five flush-mounted pressure transducers (Kulite XCQ-062-50A) are installed on the lower wall of the diffuser. In addition, five static pressure taps connected to a pressure scanner (Pressure Systems Inc., ESP-64 HD-DTC) are installed in line with the Kulite transducers. 2.2. Experimental facility and conditions The wind tunnel test was conducted using the supersonic wind tunnel (SWT) at the Chofu Aerospace Center of JAXA, which is a blowdown wind tunnel having a test section with a cross-section of 1 m by 1 m. Initially, the wind tunnel started with a constant Mach number of 2.0, and the flow plug behind the diffuser

JID:AESCTE AID:105523 /FLA

[m5G; v1.261; Prn:8/11/2019; 10:22] P.3 (1-10)

J. Yamamoto et al. / Aerospace Science and Technology ••• (••••) ••••••

3

Fig. 2. Wind tunnel test. The major dimensions in (a) are listed in Table 1. Among the five Kulite transducers, the data from the PS3 sensor mounted at x ps = 115 mm and b ps = 20 mm was used in the present research. Table 1 Main dimensions of the inlet model. The definitions of the symbols are given in Fig. 2(a). Symbol Value

zi [mm]

zo [mm]

L [mm]

bd [mm]

br [mm]

α

[deg] 10

30

40

240

60

200

70

δ

[deg]

was fully opened. Then, the flow rate through the inlet was constricted by moving the flow plug forward in a stepwise manner. The movement of the flow plug was started 2 s after the wind tunnel started. The inlet was fully closed 25 s after the wind tunnel started. Finally, the exit gradually opened again after 25 s. Schlieren images around the wedge and the entrance of the diffuser during this constriction were taken by a high-resolution still camera (Nikon D-2) with an exposure time of 1/8,000 s. The unsteady pressures measured by the five Kulite transducers were recorded at a rate of 2 kHz. The total pressure and temperature during this test were approximately p 0 = 287 kPa and T 0 = 220 K, respectively. 2.3. Experimental results A typical experimental result in this wind tunnel test is shown in Fig. 3, where Figs. 3(a) through (d) show the time histories of the horizontal displacement of the flow plug x p , the unsteady pressure measured at port PS3, the RMS, and the spectrogram of the pressure fluctuation shown in (b), respectively. The position of the PS3 sensor is shown in Fig. 2(a). The result indicates that the static pressure gradually increases with the movement of the flow plug. The RMS of the pressure fluctuation abruptly increases at 20.5 s when the displacement x p is reduced from 5 mm to 4 mm. The dominant frequency of the pressure fluctuation is 273 Hz, which is approximately constant, while substantial pressure fluctuation

Fig. 3. Experimental results. (a) Horizontal displacement of flow plug x p , (b) Pressure measured at port PS3, (c) RMS of pressure fluctuation, and (d) Spectrogram (color). Here, x p = 0 indicates that the exit area of the diffuser is fully closed. The width of the time window for obtaining (c) and (d) is 0.128 s. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

is observed. Note that this large-amplitude pressure fluctuation is preceded by a small-amplitude pressure fluctuation from approximately 18 s, the dominant frequency of which is also approximately 270 Hz. Snapshots of the schlieren images at some instants are shown in Fig. 4. The shear layer flows in the diffuser after t = 11.22 s (Fig. 4(a)). At this moment, however, no remarkable pressure fluctuation is observed, as shown in Fig. 3(c). The flow passing through the oblique and terminal shock waves are designed to reduce the total pressure to 88% of the uniform flow. As well, the relative total pressure of the flow passing through only the terminal shock wave is 72%. On the basis of the Ferri criterion, Fischer et al. [29] determined that the buzz occurs when the difference between these two relative total pressures is larger than 7%. Since the difference is 16% in our case, the Ferri (Fischer) criterion is not applicable to the onset prediction of this experiment. The shock system is fairly stable after ingesting the shear layer until t = 20.38 s (Fig. 4(b)). During this stable shock system period, the position of the shock wave triple point (STP) gradually shifts upward as the exit area of the diffuser is reduced. Subsequently, the large-amplitude oscillation of the shock system is observed (Figs. 4(c) and 4(d)) until t = 28.34 s, when the exit of the diffuser is opened again. The pressure-time histories around the snapshots (c) and (d) are displayed in Figs. 4(e) and (f). The pressure is maximum, if the final shock wave is close to the dif-

JID:AESCTE

AID:105523 /FLA

4

[m5G; v1.261; Prn:8/11/2019; 10:22] P.4 (1-10)

J. Yamamoto et al. / Aerospace Science and Technology ••• (••••) ••••••

Fig. 4. Schlieren images at some instants. (a) At the beginning of the shear layer ingestion (11.22 s), (b) just before the onset of buzz (20.38 s), and (c) and (d) during buzz (c: 23.96 s, d: 24.16 s). Pressure-time histories around t = 23.96 s and t = 24.16 s are displayed. The red broken lines are the times of snapshot (c) and (d).

fuser inlet. The pressure is minimal when the final shock wave is far from the inlet. 3. Numerical analysis 3.1. Model used for computation A modified inlet model shown in Fig. 5 was used in the numerical analysis. All the stays for installation in the wind tunnel were removed from the wind tunnel model. In order to compute the flow stably, the shape of the wedge was changed from a hollow pentagon (see Fig. 2(b)) to a filled hexagon, and the ends of diffuser walls were rounded. These modifications did not inherently affect the characterization of the inlet buzz. The wedge was used as a generator of an oblique shock wave, the structure of which was independent of the aft body shape of the wedge. The inlet buzz was mainly related to the flow rate through the diffuser [3, 10], so that the end shape of the wall did not play an important role if the flow through the diffuser was appropriately simulated. Also, in the wind tunnel model, there was a narrow gap (width 0.5 mm) between the movable upper wall and the side wall of the diffuser. In order to take the bleed flow from the diffuser through the gap into account, we also set a slit in the computational model. The major dimensions are defined in Fig. 5(a). The capture area A c is equal to the area bd zi in Fig. 2, the size of which is listed in Table 1. In addition, A e is the exit area of the diffuser, which is equal to the minimum area between the upper (or lower) wall of the diffuser and the surface of the flow plug. The cross-sectional area of the slit A b is equivalent to that of the gap in the wind tunnel model. The angle of the oblique shock wave θ is a function of wedge angle δ and the Mach number M. The gap height  y

Fig. 5. Model used for computation.

is defined as the vertical distance between the lip of the lower wall (cowl lip) and the extension of the oblique shock wave. In this calculation, θ = 40◦ with M = 2.0 and δ = 10◦ , A b / A c = 0.2,

JID:AESCTE AID:105523 /FLA

[m5G; v1.261; Prn:8/11/2019; 10:22] P.5 (1-10)

J. Yamamoto et al. / Aerospace Science and Technology ••• (••••) ••••••

5

and  y / zi = 0.11, according to the placement of the model in the wind tunnel test. 3.2. Numerical method and conditions Numerical simulation of the unsteady compressible viscous flow was conducted using the FAST Aerodynamic Routines (FaSTAR) developed by JAXA [30]. We selected URANS equations with the SA-R model [31,32] as the turbulence model from various options in FaSTAR. These equations are commonly used in industrial applications in aerospace engineering. Second-order accuracy in spatial differentiation was secured using the Monotone Upwind Scheme for Conservation Laws (MUSCL) method [33], in which an original van Leer-type flux limiter was used to improve the numerical stability. We chose the Harten, Lax, van Leer, Einfeldt, and Wada (HLLEW) scheme [34] for the Riemann solver for evaluation of the inviscid flux and the G-G based weighted least square (GLSQ) method [35] to calculate the spatial gradients of variables. Unsteady computation was performed by a second-order implicit time integration method combining the LU-SGS method [36] and the dual time stepping method [37]. We computed the flow around the inlet in the half-span region under the assumption of symmetry with respect to the xy crosssection passing through the center of inlet (see Fig. 5(b)). An H-H type hexahedral structure grid is used in the discretization of the space. We calculated the flow after the impulsive start conditions, in which all of the variables in the computational domain were initially set at the values in a uniform flow.

Fig. 6. Time history of unsteady pressure of the diffuser at the same position as the PS3 Kulite transducer installed in the wind tunnel test. Two computational grids, the total numbers of which are 3 × 106 (3M) and 24 × 106 (24M), are used in this calculation.

3.3. Grid convergence In order to examine the grid convergence of this simulation, we conducted a preliminary computation with two computational grids, the total numbers of which are 3 × 106 (3M) and 24 × 106 (24M), respectively. The minimum spacing of the grid was set at 3.0 × 10−6 m for both cases. The time step was maintained constant at 1.98 × 10−7 s. Fig. 6 shows the time history of unsteady pressure measured at the surface of the lower wall when the large-amplitude inlet buzz is observed (x p = 1 mm, A e / A c = 0.039, A b / A c = 0). The measurement point of this unsteady pressure corresponds to the position of the PS3 Kulite transducer, the position of which is indicated in Fig. 2(a). Fig. 6 reveals that the two profiles of time history agree well with each other, so that the numerical results with the 3M grids are sufficiently reliable in this simulation. Hereinafter, all of the numerical results are calculated using the 3M grids. Fig. 7 shows the schlieren image obtained by CFD calculation superimposed on the experimental schlieren photograph at t = 19.98 s where the position of flow plug is x p = 5 mm ( A e / A c = 0.191), the slit opening ratio A b / A c = 0.2 and the gap height  y / zi = 0.11 for the calculation. The figure clearly indicates that the shock system just before the onset of large-amplitude buzz is accurately simulated by our CFD calculation. 4. Results and discussion 4.1. Cycle of buzz We examine whether the numerical simulation can capture the large-amplitude cycle of inlet buzz by comparing the numerical results with the wind tunnel test data. Fig. 8 shows the time history of unsteady pressure, where the maximum amplitude of the pressure fluctuation is observed. The measurement point is the same as in Fig. 3. In the numerical simulation, the parameters are A e / A c = 0.039, A b / A c = 0.2, and  y / zi = 0.11. In the wind tunnel test, the pressure recorded 24.5 s

Fig. 7. Comparison of schlieren images of the shock system obtained by CFD (color) and wind tunnel experiment (gray scale). The position of flow plug is x p = 5 mm for both cases.

Fig. 8. Time histories of pressure by numerical simulation and wind tunnel test with x p = 1 mm ( A e / A c = 0.039).

after the flow plug is started at x p = 1 mm. The exit area is the same as for the numerical simulation ( A e / A c = 0.039). The time history shown in Fig. 8 indicates that the numerical simulation agrees well with the wind tunnel test data. The mean and the dominant frequency of the pressure fluctuation are listed in Table 2. The differences between the simulation and the wind tunnel test in terms of the mean pressure as well as the dominant frequency are acceptably small to proceed to the quantitative prediction of the inlet buzz. Fig. 9 shows pseudo schlieren profiles at some instants around the inlet during a cycle of buzz. The calculating conditions are

JID:AESCTE

AID:105523 /FLA

[m5G; v1.261; Prn:8/11/2019; 10:22] P.6 (1-10)

J. Yamamoto et al. / Aerospace Science and Technology ••• (••••) ••••••

6

Fig. 9. Pseudo schlieren profiles at some instants around the inlet during a cycle of buzz for A e / A c = 0.039. The times are assumed to be (a) 0 ms, (b) 0.77 ms, (c) 1.02 ms, (d) 2.08 ms, (e) 2.97 ms, and (f) 3.47 ms.

Fig. 10. Variety of pressure fluctuation due to difference in A e / A c . Table 2 Characteristic quantities of inlet buzz obtained by numerical simulation using a wind tunnel test. Parameter

Simulation

Experiment

Error %

Mean Pressure [kPa] Frequency [Hz]

167 276

166 279

0.5 1

the same as the result shown in Fig. 8. Fig. 9(a) shows the profile at the moment when the detached shock wave is on the most upstream side. Then, as shown in Fig. 9(b), a compression wave propagates downstream in the diffuser as the detached shock wave moves downstream. This compression wave is reflected at the exit of the diffuser, like a closed end reflection (Fig. 9(c)). The detached shock wave continuously moves downstream while the compression wave reciprocates in the diffuser. The compression wave arriving at the entrance of the diffuser interacts with the detached shock wave and is then reflected as an open end reflection (Fig. 9(d)). As a result, the detached shock wave is pushed away from the entrance of the diffuser and a expansion wave propagates downstream in the diffuser (Fig. 9(e)). Finally, the detached

shock wave moves upstream (Fig. 9(f)) until it reaches the most upstream position. Fig. 10 shows the time histories of unsteady pressure in the diffuser for four exit area A e conditions. The measurement point of these histories is the same as the position of the PS3 sensor, which is the same as in Fig. 8. Fig. 10(a) shows that the pressure is stable for the largest area condition A e / A c = 0.191, which corresponds to the horizontal position of the flow plug, i.e., x p = 5 mm. The pressure fluctuation occurs when the exit area is reduced to A e / A c = 0.153 (x p = 4 mm) (Fig. 10(b)). Then, the fluctuation increases as the exit area decreases (Figs. 10(c) and 10(d)). 4.2. Onset of buzz We characterize the inlet buzz using the shock wave triple point (STP) shown in Fig. 11, where the oblique shock wave from the wedge encounters the detached shock wave generated in front of the diffuser. The use of the STP is suitable for comparison between the wind tunnel test and the numerical simulation, because its position is clearly identified both in the wind tunnel test and in the numerical simulation. The vertical position ( y) of the STP is closely related to the flow rate through the diffuser. The flow rate

JID:AESCTE AID:105523 /FLA

[m5G; v1.261; Prn:8/11/2019; 10:22] P.7 (1-10)

J. Yamamoto et al. / Aerospace Science and Technology ••• (••••) ••••••

7

Fig. 11. Definition of STP.

through the diffuser is maximum when the position of the STP is close to the lip of lower wall of the diffuser. The vertical position of the STP moves upward as the flow rate through the diffuser decreases. For the purpose of comparison, we measured the vertical position ( y) of the STP from the wind tunnel test and the numerical simulation. The position of the STP was determined by the digital image processing technique described in Appendix. As shown in Fig. 11, the vertical position of the STP being equal to the height of the lip of the lower wall of the diffuser is indicated as 0%, and the vertical position of the STP being equal to the height of the lip of the upper wall of the diffuser is indicated as 100%. Fig. 12 shows the relation between the vertical position of the STP and the root-mean-square (RMS) of pressure fluctuation measured at the PS3 sensor whose time history has already been displayed in Fig. 3(c). The symbols represent the status of the shock oscillation: a square denotes without oscillation, a circle denotes fully developed buzz, and a triangle denotes oscillation with a relatively small amplitude. Note that the positions of the STP obtained from the wind tunnel test (open symbols) are instantaneous (not time-averaged) data taken from the time-series images of the video camera, while the positions obtained from the numerical simulation (filled symbols) are time averaged. The RMS of pressure fluctuation is close to zero if the vertical position of the STP is below 40%. The RMS of pressure fluctuation increases substantially when the vertical position of the STP is above 45%. This means that the onset of buzz is successfully determined using the vertical position of the STP. The critical positions of the STP for the onset of buzz (∼ 45%) agree well with each other in the wind tunnel test and the numerical simulation. 5. Modeling of inlet buzz using a delay differential equation The results in the preceding section strongly suggest that the onset of inlet buzz is closely related to the mass flow rate through the diffuser. Based on this perspective, we derived a simple mathematical model for buzz with the data obtained by the numerical simulation. 5.1. Mathematical model In order to simplify the situation, we considered the balance of mass in the diffuser. Let us set a control volume (CV) covering the rear part of the diffuser, as shown in Fig. 13. Based on the conservation law for mass, the mass in the CV, m(t ), satisfies the following:

dm(t ) dt

˙ in (t ) − m ˙ out (t ), =m

(1)

˙ in and m ˙ out are the mass flow rates flowing in the CV from where m the upstream and flowing out to the downstream, respectively. ˙ out consists of the mass flux through the rear The outflow rate m ˙ e and the flux due to the bleed m ˙ b . Since the end of the diffuser m Mach number of the flow through the CV is less than 0.2, the flow

Fig. 12. STP obtained experimentally and by numerical simulation.

Fig. 13. Control volume and mass flux in the diffuser.

is incompressible so that the outflow rate is linearly proportional to the pressure as well as the mass in the CV. This implies that the ˙ out should be synchronized in phase. time variations of m and m ˙ in is expected to decrease with the increase The inflow flux m in mass in the CV, because the change in the mass in the diffuser mainly comes from the excessive inflow mass flux during the cycle of buzz. However, as shown in Fig. 14, there is a phase ˙ in . The steep difference between the time variations of m and m ˙ in at t = 0.75 ms is caused by arrival of the compresincrease in m sion wave propagating downstream of the diffuser (see Fig. 9(b)). The increase of mass in CV occurs late to the increase of mass ˙ in from t = 1.9 ms to t = 3.4 ms is flow rate. The decrease in m caused by the expansion wave (see Fig. 9(d)–(f)). The decrease in the mass is behind the decrease in flow rate. A cross-correlation analysis indicates that the phase of m is approximately 1.2 ms be˙ in if we assume that their changes are out of hind the phase of m phase with each other. We denote this time delay by τ . ˙ out (t ) and beFig. 15 shows Lissajous plots between m(t ) and m ˙ in (t + τ ). The outflow rate m ˙ out (t ) increases as tween m(t ) and m the mass m increases. As shown by the fitted lines in Fig. 15, the ˙ out is linear at a constant exit area A e : relation between m and m

˙ out (t ) = C out (m(t ) − mn ) , m

(2)

where mn is a constant, and the slope of fitted lines C out decreases as the exit area decreases. This is because the mass flow rate through the flow plug is in proportion to the mass due to its ˙ in is approxchoked flow [38]. The relationship between m and m imately linear as expected time histories of these variables shown in Fig. 14, albeit with some variation in the data. The negative ˙ in in Fig. 15 is also reasonable because the increase in slope of m mass (as well as pressure) in the diffuser apparently suppresses to

JID:AESCTE

AID:105523 /FLA

[m5G; v1.261; Prn:8/11/2019; 10:22] P.8 (1-10)

J. Yamamoto et al. / Aerospace Science and Technology ••• (••••) ••••••

8

with the largest C out is attenuated, the result with the second largest C out is neutral, and the other two results are amplified. The dominant frequency of these solutions is 283 Hz, which is fairly close to the frequency of the buzz cycle obtained from the experiment and the numerical simulation (see Table 2). In conclusion, the delay differential equation modeled in the present study (Eq. (4)) describes not only the bifurcation of the oscillatory solutions but also the frequency of the buzz in quantitative manner. A minor problem is that there is no limit cycle in the solution of Eq. (4) because we assume that the system is linear. In order for the solution of an equation to include a limit cycle, it is necessary to include a non-linear term in the equation [40]. ˙ in (red line) and m (blue line), where A e = 3.8% (x p = Fig. 14. Time histories of m 1 mm) whose contour of the density gradient has been displayed in Fig. 9. Points (a)-(f) on the red line indicate the time of the pseudo schlieren images in Fig. 9. The broken line indicates the mass time history delayed by 1.2 ms.

5.2. Criteria for the onset of buzz In order to determine the stability of the solution of Eq. (4), we seek (complex) characteristic root λ such that m(t ) = exp(λt ) is a solution of Eq. (4). The characteristic equation is

λ = C in exp (−τ λ) − C out .

(5)

The complex root λ = α + i β (α : growth rate, β : angular frequency) is substituted into Eq. (5), and then the equation is separated into real and imaginary parts using exp(−i θ) = cos θ − i sin θ as

α = C in exp (−τ α ) cos (τ β) − C out ,

(6)

β = −C in exp (−τ α ) sin (τ β) .

(7)

We focus on finding the critical point of the system that is unstable, namely having a growth rate of α = 0. In this case,

0 = C in cos (τ βc ) − C out ,c ,

βc = −C in sin (τ βc ) ,

˙ out (t ) and between (b) m(t ) and Fig. 15. Lissajous plots between (a) m(t ) and m ˙ in (t + τ ). m

ingest the flow from the entrance of diffuser. These data approximately collapse into a line, irrespective of the exit area:

˙ in (t ) = C in (m(t − τ ) − m0 ) , m

(3)

where C in and m0 are constant. The fitted values for the parameters C in , C out , m0 , and mn are listed in Table 3. Substituting Eqs. (2) and (3) into Eq. (1), we obtain a delay differential equation in terms of m:

dm(t ) dt

= C in (m(t − τ ) − m0 ) − C out (m(t ) − mn ) .

(4)

The delay differential equation has a wide variety of solutions [39]. For example, if the delay τ is zero, Eq. (4) results in a first-order ordinary differential equation, which does not have an oscillatory solution. The system in our case contains a delayed negative feedback due to C in < 0 and τ > 0, which may induce an oscillatory solution. The inlet buzz is a phenomenon that is peculiar to a compressible flow, because the time delay is caused by the finite speed of the compression wave and the expansion wave in the diffuser. Fig. 16 shows the numerical solution of Eq. (4) using the parameters listed in Table 3. These solutions are obtained by the fourth-order Runge-Kutta method. The mass at t = 0 ms is determined using the density of the uniform flow. The graphs indicate that all the solutions exhibit oscillation. Among them, the result

(8) (9)

where βc and C out ,c are the critical values for β and C out , respectively. Substituting τ = 1.20 ms and C in = −2.14 × 103 s−1 into the above equations, we obtain

βc = 1.79 × 103 ,

C out ,c = 1.12 × 103 .

(10)

The value of C out ,c is appropriate from the calculation results shown in Fig. 16, where the solution for C out = 1.16 × 103 s−1 exhibits the oscillation of constant amplitude. The critical frequency f c (= βc /(2π )), which corresponds to the frequency of buzz, is 285 Hz. This value is also fairly close to the frequency of the solutions shown in Fig. 16. We can evaluate the onset and dominant frequency of supersonic inlet buzz if we obtain the sensitivities C in and C out and the time delay τ from the numerical simulation or wind tunnel test. 6. Conclusions We conducted a comprehensive study about the onset of supersonic inlet buzz. A numerical simulation using a standard URANS approach was compared with a wind tunnel test for an externalcompression inlet model. Through the comparison, we show that our numerical simulation is sufficiently accurate for the prediction of the inlet buzz flow. An analysis of the position of the STP indicates that the onset of inlet buzz is closely related to the mass flow rate through the diffuser. Based on this perspective, we derived a simple delay differential equation (Eq. (4)) where we consider the balance of mass at the rear part of the diffuser. The derived equation contains a delayed negative feedback caused by the presence of the compression and expansion waves propagating in the diffuser with a finite

JID:AESCTE AID:105523 /FLA

[m5G; v1.261; Prn:8/11/2019; 10:22] P.9 (1-10)

J. Yamamoto et al. / Aerospace Science and Technology ••• (••••) ••••••

9

Table 3 Summary of the parameters for the buzz model (Eq. (4)). Status

x p [mm]

C out [×103 s−1 ]

mn [×10−3 kg]

C in [×103 s−1 ]

m0 [×10−3 kg]

τ [ms]

buzz buzz buzz transition stable

1 2 3 4 5

0.746 0.876 1.08 1.16 1.26

2.77 × 10−2

−2.14

0.198

1.20

Fig. 16. Numerical solution of the buzz model (Eq. (4)) by the 4th-order Runge-Kutta method. The dominant frequency of these solutions is 283 Hz.

speed. The presence of this feedback mechanism is the origin of the self-sustained oscillation of the buzz. The equation is characterized by three parameters: the delay time τ due to the propagation of pressure waves in the diffuser, the sensitivity of the inflow rate to the mass C in , and the sensitivity of the outflow rate to the mass C out . A calculation using the parameters determined from the numerical simulation indicates that the onset and the dominant frequency of the buzz are predicted quantitatively well by this delay differential equation. The critical values of the parameters for the onset are derived by the stability analysis for the equation. Declaration of competing interest

Fig. 17. Image processing procedure for determination of the STP. Dashed lines indicate the moving range of the shock wave.

There is no competing interest for authors in this paper.

We received a grant for collaborative research from the JAXA Aviation Program Group (2010–2012) as a startup budget. We would like to thank our graduate students, S. Kimura, H. Hori, Y. Tsuchimoto, Y. Saito, M. Doi who conducted the analysis of the wind tunnel test data and the numerical simulation. We used JAXA Supercomputer System generation 2 (JSS2) for the numerical simulation.

schlieren images captured by the high-speed camera in the wind tunnel test and the pseudo schlieren images obtained by the numerical simulation. A schematic diagram of the technique is shown in Fig. 17. The detached shock wave was highlighted by binarization of the original image using open image processing software (Image J [41]). Then, we calculated the center of mass of the highlighted pixels along the line of oblique shock wave within a thin bounding box covering its entire movement range. The vertical position ( y) of this center of mass was defined as the position of the STP in this image.

Appendix. Determination of the shock wave triple point

References

Acknowledgements

The vertical position ( y) of the STP both from the wind tunnel test and the numerical simulation was measured by a digital image processing technique. The technique was applied to the

[1] M. Honda, K. Yoshida, D-SEND # 2 flight demonstration for low sonic boom design technology, in: Proceedings of the 29th Congress of the International Council of the Aeronautical Sciences, 2014, 2014-0505.

JID:AESCTE

10

AID:105523 /FLA

[m5G; v1.261; Prn:8/11/2019; 10:22] P.10 (1-10)

J. Yamamoto et al. / Aerospace Science and Technology ••• (••••) ••••••

[2] P. Coen, A. Loubeau, B. Pauer, NASA’s low boom flight demonstration: assessing community response to supersonic overflight of quiet supersonic aircraft, J. Acoust. Soc. Am. 141 (5) (2017) 3624, https://doi.org/10.1121/1.4987783. [3] J. Seddon, E.L. Goldsmith, Intake Aerodynamics, 2nd ed., AIAA, 1999. [4] J.W. Slater, SUPIN: A computational tool for supersonic inlet design, in: 54th AIAA Aerospace Sciences Meeting, San Diego, California, USA, 4–8 January 2016, AIAA 2016-0530, 2016. [5] L.M. Berra, J.W. Slater, S.M. Olcmen, Conceptual redesign of the B-1B bomber inlets for improved supersonic performance, Aerosp. Sci. Technol. 45 (2015) 476–483, https://doi.org/10.1016/j.ast.2015.06.017. [6] T.R. Conners, D.C. Howe, Supersonic inlet shaping for dramatic reductions in drag and sonic boom strength, in: 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, USA, 9–12 January 2016, AIAA 2006-0030, 2006. [7] A. Murzionak, J. Etele, Rapid supersonic performance estimation for a novel RBCC engine inlet, Aerosp. Sci. Technol. 32 (1) (2014) 51–59, https://doi.org/10. 1016/j.ast.2013.12.007. [8] E.L. Goldsmith, J. Seddon, Practical Intake Aerodynamic Design, AIAA, 1993. [9] R.W. Newsome, Numerical simulation of near-critical and unsteady, subcritical inlet flow, AIAA J. 22 (10) (1984) 1375–1379, https://doi.org/10.2514/3.48577. [10] P.-J. Lu, L.-T. Jain, Numerical investigation of inlet buzz flow, J. Propuls. Power 14 (1) (1998) 90–100, https://doi.org/10.2514/2.5254. [11] C.L. Dailey, Supersonic diffuser instability, J. Aeronaut. Sci. 22 (11) (1955) 733–749, https://doi.org/10.2514/8.3452. [12] S. Trapier, S. Deck, P. Duveau, Delayed detached-eddy simulation and analysis of supersonic inlet buzz, AIAA J. 46 (1) (2008) 118–131, https://doi.org/10. 2514/1.32187. [13] S. Trapier, P. Duveau, S. Deck, Experimental study of supersonic inlet buzz, AIAA J. 44 (10) (2006) 2354–2365, https://doi.org/10.2514/1.20451. [14] J.Y. Oh, F. Ma, S.-Y. Hsieh, V. Yang, Interactions between shock and acoustic waves in a supersonic inlet diffuser, J. Propuls. Power 21 (3) (2005) 486–495, https://doi.org/10.2514/1.9671. [15] X. Shi, J.-T. Chang, W. Bao, D. Yu, B. Li, Supersonic inlet buzz margin control of ducted rockets, Proc. Inst. Mech. Eng., Part G, J. Aerosp. Eng. 224 (10) (2010) 1131–1139, https://doi.org/10.1243/09544100JAERO687. [16] T. Herges, J. Dutton, G. Elliott, High-speed schlieren analysis of buzz in a relaxed-compression supersonic inlet, in: 48th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit Atlanta, Georgia, USA, 30 July 2012 - 01 August 2012, AIAA 2012-4146, 2012. [17] E. Kwak, S. Lee, Numerical study of the effect of exit configurations on supersonic inlet buzz, in: 31st AIAA Applied Aerodynamics Conference, San Diego, CA, USA, June 24-27, 2013, AIAA 2013-3025, 2013. [18] W. Hong, C. Kim, Computational study on hysteretic inlet buzz characteristics under varying mass flow conditions, AIAA J. 52 (7) (2014) 1357–1373, https:// doi.org/10.2514/1.J052481. [19] S. Candon, E. Loth, M. Rybalko, S. Hirt, Acoustically induced shock oscillations in a low-boom inlet, AIAA J. 54 (7) (2016) 2134–2148, https://doi.org/10.2514/ 1.J054300. [20] M.R. Soltani, J. Sepahi-Younsi, Buzz cycle description in an axisymmetric mixed-compression air intake, AIAA J. 54 (3) (2016) 1040–1053, https://doi. org/10.2514/1.J054215. [21] H. Chen, H.-J. Tan, Q.-F. Zhang, Y. Zhang, Buzz flows in an external-compression inlet with partially isentropic compression, AIAA J. 55 (12) (2017) 4286–4295, https://doi.org/10.2514/1.J056066. [22] H. Chen, H.-J. Tan, Q.-F. Zhang, Y. Zhang, Throttling process and buzz mechanism of a supersonic inlet at overspeed mode, AIAA J. 56 (5) (2018) 1953–1964, https://doi.org/10.2514/1.J056674.

[23] W. Shi, J. Chang, Y. Wang, W. Bao, X. Liu, Buzz evolution process investigation of a two-ramp inlet with translating cowl, Aerosp. Sci. Technol. 84 (2019) 712–723, https://doi.org/10.1016/j.ast.2018.11.016. [24] J. Chang, D. Yu, W. Bao, C. Wang, T. Chen, Mathematical modeling and rapid recognition of hypersonic inlet buzz, Aerosp. Sci. Technol. 23 (1) (2012) 172–178, https://doi.org/10.1016/j.ast.2011.07.008. [25] M. Farahani, A. Daliri, J.S. Younsi, Supersonic inlet buzz detection using pressure measurement on wind tunnel wall, Aerosp. Sci. Technol. 86 (2019) 782–793, https://doi.org/10.1016/j.ast.2019.02.002. [26] M.R. Soltani, M. Farahani, Effects of angle of attack on inlet buzz, J. Propuls. Power 28 (4) (2012) 747–757, https://doi.org/10.2514/1.B34209. [27] M.R. Soltani, A. Daliri, J.S. Younsi, M. Farahani, Effects of bleed position on the stability of a supersonic inlet, J. Propuls. Power 32 (5) (2016) 1153–1166, https://doi.org/10.2514/1.B36162. [28] A. Ferri, L.M. Nucci, The Origin of Aerodynamic Instability of Supersonic Inlet at Subcritical Conditions, NACA Research Memorandum L50K30, 1951. [29] S.A. Fisher, M.C. Neale, A.J. Brooks, On the sub-critical stability of variable ramp intakes at Mach numbers around 2, in: Aeronautical Research Council, Research and Memoranda No. 3711, 1970. [30] A. Hashimoto, K. Murakami, T. Aoyama, K. Ishiko, M. Hishida, M. Sakashita, P. Lahur, Toward the fastest unstructured CFD code “FaSTAR”, in: 50th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, Nashville, Tennessee, USA, 9 – 12 January 2012, AIAA 2012-1075, 2012. [31] J. Dacles-Mariani, G.G. Zilliac, J.S. Chow, P. Bradshaw, Numerical/experimental study of a wingtip vortex in the near field, AIAA J. 33 (9) (1995) 1561–1568, https://doi.org/10.2514/3.12826. [32] Z. Lei, Effect of RANS turbulence models on computation of vortical flow over wing-body configuration, Trans. Jpn. Soc. Aeronaut. Space Sci. 48 (161) (2005) 152–160, https://doi.org/10.2322/tjsass.48.152. [33] B. van Leer, Towards the ultimate conservative difference scheme. V. A secondorder sequel to Godunov’s method, J. Comput. Phys. 32 (1) (1979) 101–136, https://doi.org/10.1016/0021-9991(79)90145-1. [34] S. Obayashi, Y. Wada, Practical formulation of a positively conservative scheme, AIAA J. 32 (5) (1994) 1093–1095, https://doi.org/10.2514/3.12104. [35] E. Shima, K. Kitamura, T. Haga, Green-Gauss/weighted-least-squares hybrid gradient reconstruction for arbitrary polyhedra unstructured grids, AIAA J. 51 (11) (2013) 2740–2747, https://doi.org/10.2514/1.J052095. [36] S. Yoon, A. Jameson, Lower-upper symmetric-Gauss-Seidel method for the Euler and Navier-Stokes equations, AIAA J. 26 (9) (1998) 1025–1026, https://doi.org/ 10.2514/3.10007. [37] C. Hirsch, Numerical Computation of Internal and External Flows, 2nd edition, Butterworth-Heinemann, Oxford, 2007. [38] J.D. Anderson, Modern Compressible Flow: With Historical Perspective, 3rd edition, McGraw Hill, New York, 2003. [39] H. Smith, An Introduction to Delay Differential Equations with Applications to the Life Sciences, Springer-Verlag, New York, 2011. [40] S.H. Strogatz, Nonlinear Dynamics and Chaos with Applications to Physics, Biology, Chemistry, and Engineering, second edition, CRC Press, Boca Raton, 2015. [41] C.A. Schneider, W.S. Rasband, K.W. Eliceiri, NIH image to ImageJ: 25 years of image analysis, Nat. Methods 9 (2012) 671–675, https://doi.org/10.1038/nmeth. 2089.