Large-eddy simulations of flow and aeroacoustics of twin square jets including turbulence tripping

In this study, we investigate the flow and aeroacoustics of twin square (i.e., aspect ratio of 1.0) jets by implicit large-eddy simulations (LESs) under a nozzle pressure ratio of 3.0 and a temperature ratio of 1.0 conditions. A second-order central scheme coupled with a modified Jameson ’ s artificial dissipation is used to resolve acoustics as well as to capture discontinuous solutions, e.g., shock waves. The flow boundary layer inside of the nozzle is tripped, using a small step in the convergent section of the nozzle. The time-averaged axial velocity and turbulent kinetic energy of LES with boundary layer tripping approaches better to particle image velocimetry experimental data than the LES without turbulence tripping case. A two-point space – time cross-correlation analysis suggests that the twin jets are screeching and are coupled to each other in a symmetrical flapping mode. Intense pressure fluctuations and standing waves are observed between the jets. Spectral proper orthogonal decomposition (SPOD) confirms the determined mode and the relevant wave propagation. The upstream propagating mode associated with the shock-cell structures is confined inside jets. Far-field noise obtained by solving Ffowcs Williams and Hawkings equation is in good agreement with the measured acoustic data. The symmetrical flapping mode of twin jets yields different levels of the screech tone depending on observation planes. The tonalities — the fundamental tone, second and third harmonics — appear clearly in the far-field, showing different contributions at angles corresponding to the directivities revealed by SPOD.


I. INTRODUCTION
For past several decades, supersonic jets have been used as an essential way to generate thrust of high-speed flight vehicles.The interaction of the shock-cell structure-which is generated by the pressure difference between supersonic jets and the ambient air-jet shear layers, and turbulent flows induces surrounding pressure fluctuations that generate significantly intense acoustics at near-and far-fields.The large-scale turbulent mixing noise and broadband shock-associated noise (BBSAN) contribute to the broadband noise spectra.On the other hand, the screech noise, which can be measured upstream of jet flows, is generated by the closed feedback loop of coherent flow structures convecting downstream in jet shear layers and upstream propagating acoustic waves induced by perturbed shock cells.
As a mean to achieve a higher performance thrust power and/or a precautionary reason for one jet engine failure, multiple jet configurations have been widely adopted in modern aircrafts.However, coupling modes of twin jets, yielding flow and acoustic resonances, can appear when two jets are operated in proximity to each other.In particular, it was experimentally found that a symmetrical (with respect to the symmetry plane of the twin nozzle) coupling mode causes significantly amplified levels of the screech noise. 1,24][5][6][7][8] In addition, passive 1,3 and active 7 flow control methods for decoupling twin jets were also assessed.By using a two-dimensional vortex sheet model, Tam and Seiner 9 early suggested that the symmetrical and antisymmetrical flapping modes are the possible coupling modes of the twin circular jet configuration.Recently, Knast et al. 10 experimentally visualized more types of coupling modes including the flapping and helical ones by applying a two-point cross-correlation analysis.
Due to the requirement of huge computations and sensitivity for capturing the screech phenomenon of coupled twin jets, numerical simulations dedicated to the twin jets were initiated recently.The development of twin jets and the jet-to-jet interactional influences were investigated by Goparaju and Gaitonde 11 using large-eddy simulation (LES).Including a far-field noise prediction by solving Ffowcs Williams and Hawkings equation, Erwin et al. 12 performed LES for the hot twin jets impinging on a deflector.Gao et al. 13 first simulated the symmetrical flapping mode of twin circular jets.The flows and near-field noise spectra of the twin jets were compared to those of the single jet case, and the frequency change due to shifting of the screech sources was assessed by applying dynamic mode decomposition (DMD).Ahn et al. 14 considered the two twin nozzle configurations having different inter-nozzle distance (i.e., nozzle-to-nozzle centerline spacing) and investigated the characteristics of near-field pressure fluctuations of the opposite flapping modes (i.e., symmetrical and antisymmetrical flapping modes) by LES.
Efforts to identify the flow and noise characteristics of a rectangular jet, which has the advantages of the airframe integration and thrust-vectoring, have been actively conducted analytically, 15,16 experimentally, [17][18][19][20][21][22][23][24] and numerically [25][26][27][28][29][30][31][32] for various aspect ratios (ARs) of the nozzles.For a twin rectangular nozzle configuration, the flow-acoustic resonance of coupled twin jets was early found, 3,33 and the dynamics of the jets represented by coupling modes and associated near-field acoustics had been examined experimentally. 34In the experimental study on the coupled twin rectangular jets of high AR of 4.0 by Raman and Taghavi 35 and Raman, 36 the flapping mode perpendicular to the twin-jet plane was investigated in detail.By quantifying the coherence between the twin jets, the switch of the symmetric and antisymmetric modes was found, which is relevant to the location of effective sources.Recently, Karnam et al. 37,38 thoroughly investigated the nature of flow and noise characteristics of coupled twin rectangular jets of AR of 2.0 by adopting signal analyses of the cross-correlation and the wavelet analyses, and image processing using spectral proper orthogonal decomposition (SPOD).For the same nozzle configuration of Karnam et al., 37,38 the development of twin rectangular jets has been estimated by Reynolds-averaged Navier-Stokes (RANS) and the result was compared to the single rectangular jet. 39The influence of wall models on the far-field noise directivity and noise spectra was also quantified by Viswanath et al.; 40 however, details of the tonal noise generated by the twin-jet coupling mode were excluded from the discussions.Most recently, Jeun et al. 41,42 successfully simulated overexpanded screeching twin rectangular jets representing the antisymmetrical flapping mode with respect to the symmetry plane of the twin nozzle, of which prediction was reliable through validations by the experimental data and SPOD analysis.
Simulations of twin-jet configurations using LES are possible today thanks to the increased computing power capabilities and computing resources.This allows us to enhance the physical understandings of the jet-to-jet interaction.The effort to explore twin rectangular jets by LES has been launched, with results associated with coupled twin rectangular jets exhausting from nozzles of AR of 2.0 just recently being published. 41,42There is a need to further investigate different rectangular nozzle configurations and the corresponding aeroacoustics characteristics.On the other hand, in experiments, disturbed jet shear layers can arise due to a turbulent boundary layer that develops inside the nozzle.Therefore, a numerical study that includes such physics and examines their influences on the generation of the resonance mode in twin jets is also necessary.
In this study, we aim to examine the aerodynamics and aeroacoustics of twin rectangular jets with AR of 1.0 (i.e., twin square jets) by LES-a nozzle pressure ratio (NPR) of 3.0 is selected for our study as we are particularly interested in investigating the flow physics and aeroacoustics under the over-expanded screeching jet condition, which has a potential to exhibit a coupling mode (i.e., resonance) for a twinnozzle configuration.The temperature ratio (TR) is fixed at 1.0 for a basic study.University of Cincinnati designed the twin square nozzle and recently performed the experiment. 43Furthermore, our preliminary LES study, which assessed the possibility of capturing the coupling mode of twin square jets, 44 has been improved to include the tripping of the boundary layer inside the twin square nozzle for better prediction of the initial (disturbed) jet shear layers and associated flow/acoustics.To the best of our knowledge, this has not been done for coupled twin jets by other researchers.In addition, data processing with SPOD is applied for revealing the dominant flow structures and pressure modes associated with the coupling phenomenon.The newly performed validations and analysis on far-field noise are also included.This paper is organized as follows: In Sec.II, details of the nozzle geometry and jet operating conditions are introduced.In Sec.III, the computational setup including the flow solver, the grid generation, and the boundary conditions is introduced.The validation of the present computational data (flow and acoustics) and the analysis on the flow and noise characteristics of coupled twin square jets are described in Sec.IV.Summary and conclusions are presented in Sec.V.

II. NOZZLE GEOMETRY AND JET OPERATING CONDITION
In this study, jets exhausting from a twin square nozzle are investigated by large-eddy simulation (LES) calculations.The global shape of the nozzle can be found in Fig. 1(a).The outer surface of the nozzle initiates in a large circular shape and gradually changes to in cross section along the streamwise direction to a smaller rectangular shape.This continues to two square shapes until the nozzle exit plane is reached.Because of this particularity, a cavity region is formed between the two (twin) nozzles.On the other hand, detail scales of the essential part of the nozzle inner surface are presented in Figs.1(b) and 1(c), which are normalized by the nozzle exit height (h) of 16.61 mm.The air-flow is supplied from the upstream pipe whose diameter is 6.61h, and is then accelerated to moderate speed until it reaches to the straight conduits.The same mass-flow rate can be expected in each conduit under an ideal scenario.Near the end of the nozzle, the air-flow is additionally accelerated to supersonic speed through a converging-diverging (C-D) section that has an area ratio of 1.168 (of the nozzle throat and the exit).It corresponds to a design Mach number of approximately 1.5 (more precisely, 1.488).It is important to note that the C-D geometry is only introduced at the upper and bottom inner walls.In other words, the side (lateral) inner walls are flat.In addition, edges at the nozzle throat are sharp; thus, generation of a set of oblique shock waves is expected here.For convenience of discussion, three reference planes are defined as the flat-side, C-D side, and symmetry planes, as presented in Fig. 1(a).
In the present LES, an over-expended jet operating condition with a nozzle pressure ratio (NPR) of 3.0 is selected along with a temperature ratio (TR) of 1.0.Details of the major parameters for the simulations are listed in Table I.Note that a Reynolds number is calculated based on a nozzle equivalent diameter (D eq ) that can be defined as 1. 25 satisfying the same pressure drop between the round and square nozzles. 45

III. COMPUTATIONAL METHODOLOGIES
In this work, LESs are performed by using an unsteady compressible finite-volume Navier-Stokes solver on an unstructured grid-based framework. 46A second-order central difference scheme is used for a spatial discretization of the non-linear convection terms, which is effective way to resolve acoustic waves by the no-dissipative feature.For the temporal integration, an explicit standard four-stage Runge-Kutta method is applied.
On the other hand, supersonic jets can generate strong shock waves.It is thus necessary to impose proper amount of numerical dissipation on local regions involving shock waves.In the present LES, a modified version 29,47 of Jameson's artificial dissipation 48 is adopted at the end of every first stage of the Runge-Kutta method.The numerical dissipation (d) applied between nodes 0 and 1 can be expressed as where U is the conservative variable, and D is a second-order differential operator (i.e., Laplacian operator).The local spectral radius (r 01 ) is defined as where u 01 is the averaged cell face speeds, c 01 is the averaged cell face speed of sound, and S 01 is the surface of the cell face.The function u 01 , which is associated with the stretching in the grid, is defined as where u 0 ¼ r 0 =4r 01 ð Þ 0:3 .The dissipation functions,

Jet operating condition NPR TR u j M j Re
Over-expanded screeching jet 3.0 1.0 399 m/s 1.36 8:85 and respectively, where and Here, p is the pressure, and m is the number of neighbors to a corresponding node.u and x are the velocity and the vorticity, respectively.s ¼ 3ðN 0 þ N 1 Þ=N 0 N 1 is a scaling factor determined by the number of neighbors N. The performance of the artificial dissipation model for typical benchmark problems is well described in the paper by Gojon et al. 49 An implicit LES approach (or monotonically integrated large-eddy simulation, MILES 50 ) is applied, which means that no explicit sub-grid-scale (SGS) model is implemented.In other words, the physical dissipation originated by small-scale eddies is implicitly considered by taking the advantage of truncation errors.Viscosity ðl) is estimated by Sutherland's law, and a fourth-order central difference is used to calculate viscous flux term.The set of aforementioned numerical schemes has been proved for many supersonic jet aeroacoustics applications, 28,29,32,[51][52][53] with good agreements between predictions and experimental data in terms of flow fields and far-field acoustics.
A three-dimensional computational grid and its midlongitudinal cut along the flat-side plane are presented in Fig. 2. The nozzle exit plane is located at z/h ¼ 0, and the size of the domain in the axial direction is 83h of which 9h occupies the upstream region.On the other hand, the domain has a length of around 15h in the radial direction at the nozzle exit plane, and the size gradually expands along the upstream and downstream directions by considering expansion of jet plumes and propagation of acoustic waves.The two inlet boundaries are set to a total pressure of 3.0 P atm and total temperature of 1.0 T atm , with zero turbulent intensity, representing a laminar flow state at the inlets.The inner and outer surfaces of the nozzle are treated with an adiabatic no-slip wall condition using a weak-formulation. 54A characteristics boundary condition using invariants is applied to the surrounding far-field boundaries, with the implementation of the ambient air state (P atm and T atm ).
In the region other than the inside of the nozzle, the smallest grid of Dx, Dy (i.e., the size perpendicular to the streamwise direction) of h/320 is used in the region where jet shear layers are expected to appear (i.e., h/2 from the jet center).The grid size is gradually increased in the streamwise and radial directions with a growth ratio controlled below 2% to prevent numerical oscillations.As a result, the size of Dz ¼ h/100, h/28, h/14, and h/7 is generated at z/h ¼ 0, 4, 15, and 40, respectively.On the other hand, a large growth ratio of nearly 5% is adopted in the sponge zone region of z/h > 40 to ensure that the flows and acoustics can be dissipated smoothly.240 grid points are populated in the azimuthal direction with respect to each square nozzle; thus, the computational grid size of around h/60 is assigned at the four edges of each square nozzle exit.
The twin square nozzle tested in the experimental facility at University of Cincinnati was made by using a three-dimensional (3D) printer.It is expected that a boundary layer will develop inside the nozzle.A few studies addressed the boundary layer development for single jet simulations by implementing an artificial geometry 55,56 or a wall model with surface roughness. 40In the current study, we applied a geometrical tripping method to our twin square nozzle configuration.A detail of the step used and of the computational mesh in the region of the step is shown in Fig. 3.The step is located in the middle of the converging section instead of the diverging section to avoid generation of shock waves at the step location due to the artificial geometry.Langenais et al. 56 implemented the step of 0.01D (D is a diameter of the nozzle exit) in their LES of a hot supersonic single round jet, and achieved turbulent intensity of around 0.04 at the nozzle exit.Since the nozzle configuration has a steep slope of the converging section, the step of 1% of the nozzle height was enough to generate proper amplitude of turbulent intensity.The height of the step in the current setup is 0.02h justified by the lower slope in the converging section as compared with the aforementioned study.To maintain a reasonable computational cost, we used the smallest grid size of h/320 at the first layer of the wall, which corresponds to a dimensionless distance to the wall y þ of 10 to 20 inside the nozzle.While a wall-layer model can be applied 40,55,57 to solve the lack of grid points near the wall in LES, we chose to use the suggested grid size without any wall-layer model, based on the acceptable results produced by the work of Langenais et al. 56 The grid resolution is insufficient to resolve the viscous sublayer but allows us to convect forced velocity fluctuations 56 and maintain a reasonable computational cost.LES results both without and with the tripping of boundary layer are discussed in this paper.Consequently, the entire computational domain consists of 145 Â 10 6 and 150 Â 10 6 grid points, respectively.Note that we named those cases as LES-NT (no tripping) and LES-WT (with tripping).
A time step of 0.0012h/u j was employed with the unsteady simulations.The acquisition of flow and acoustic data was initiated after the jet was fully developed.The calculations of LES-NT and LES-WT were performed on 960 processors.

IV. RESULTS AND DISCUSSION A. Influence of turbulence tripping on the initial jet shear layers
The turbulence tripping may induce large unsteadiness in the flow near the walls inside the nozzle.Indeed, larger flow structures can be observed just downstream of the step as depicted in Fig. 4(a), which are convected along the convergent section without significant loss of the strength.The vorticity contours in a few cross-sectional planes are depicted in Fig. 4(b).As the flow is traveling through the nozzle throat and the divergent section, the flow structures are stretched in the streamwise direction due to acceleration, losing strength.When leaving the nozzle exit, Kelvin-Helmholtz instability enhances the mixing of the jets with ambient air so that the vorticity magnitude is again increased immediately.A similar observation has been reported in a few papers for single nozzle configurations. 55,56The vorticity distribution appears to be uniform along the walls just downstream of the step (A).It becomes much irregular and weaker at the downstream sections (B and C).The lowest level of the vorticity is observed at the nozzle exit (C) among the three different locations.
As aforementioned, the purpose of the tripping is to obtain a proper level of turbulent intensity (or velocity fluctuations) at the nozzle exit.Turbulent intensity plots of both LES-NT and LES-WT are shown in Fig. 5.The tripping of the boundary layer truly enhances velocity fluctuations, in particular close to the nozzle walls.It can be observed that the tripping contributes larger fluctuations in the streamwise direction (u rms =u j Þ than in the transverse directions (v rms =u j and w rms =u j Þ.The stretched vortical flow structures demonstrated in Fig. 4, which are due to acceleration, are likely associated with this result.Consequently, the total turbulent intensity reaches around 0.06 for LES-WT case, which is slightly higher level as compared to LES of hot supersonic jet of Langenais et al. 56 (i.e., 0.04).

B. Flow characteristics of twin square jets
1. Time-averaged jet axial velocity Figure 6 contrasts the normalized time-averaged axial velocity fields obtained by particle image velocimetry (PIV) and LES.Note that the velocity is normalized by the fully expanded jet velocity u j (¼399 m/s).When leaving the nozzle exit, jet flows generate shock-cell structures owing to weak pressure with respect to the ambient air state (i.e., over-expanded jet) so that jet flows experience accelerations and decelerations as traveling downstream.It seems that the strength of the shock-cell structures by LES is slightly over-predicted as compared to the PIV experimental result.Also, a longer length of the supersonic region that is influenced by the shock-cell structures is observed in the LES results.Nevertheless, with the LES-WT, the same region is approaching close to the PIV data, the development of the jet shear layers being enhanced by tripping the boundary layer.The PIV result shows slight asymmetric velocity distribution in the C-D-side plane with respect to the jet centerline, whereas LES results do not.Nozzle geometry defects or unknown physics of twin jet interactions may cause such result, which requires further examinations.
For detailed comparisons, the time-averaged axial velocity along the jet centerline (y/h ¼ 0) is presented in Fig. 7.A locally decelerated region can be observed just before the nozzle exit (z/h ¼ 0), which is due to the weak oblique shocks generated at the sharp nozzle throat as can be seen in the subfigures on the right-hand side of Figs.6(b) and 6(c).It can be seen that the valleys of axial velocity in LES are deeper than the experimental ones, which implies the over-prediction of the shock waves in LES as compared with experiments.Nevertheless, the overall LES trends in terms of both amplitudes and phase of the flow structures are in good agreement with the experimental data.LES WT case predicts closer to experiments solution in terms of the length of the supersonic region underlying the influence of the shock-cell structures and decay of the axial velocity.Figure 8 presents the normalized axial velocity profiles at several axial locations downstream of the nozzle exit plane in two orthogonal planes, i. e., the flat-side and the C-D-side planes.Overall, there is good agreement with the experimental data.The development of the jet shear layers is better approaching the PIV when the turbulence tripping treatment is implemented (LES-WT case).

Turbulent kinetic energy (TKE)
Turbulent kinetic energy (TKE) distributions normalized by the square of fully expanded jet velocity (i.e., TKE/u 2 j ) in the flat-side plane and in the C-D side plane as measured by the PIV experiment and as predicted by LES are contrasted in Fig. 9. Also, the corresponding TKE plots along y-lateral direction at different axial locations downstream from the nozzle exit plane are also depicted in Fig. 10.The LES data are presented for both set-ups, i.e., LES-NT and LES-WT.
In the case of the experiment, the turbulent state of boundary layers inside the nozzle and the amplitude of the instability waves can be a trigger for an early generation of fine-scale turbulent flow structures, closer to the nozzle lip.This contributes to different turbulence production levels as compared with the simulations in the initial development region of the jet shear layers, closer to the nozzle trailing edge.On the other hand, a particular pattern of TKE levels is observed at the inner jet shear layers suggesting a stronger interaction between the twin jets in the experimental setup as compared with the LES predictions.Approximately at downstream locations ranging from 5 < z/h < 10, the maximum level of TKE can be observed in the lateral jet shear layers.Larger levels of TKE within the more spread, thicker jet shear layers as compared with the LES, are also shown in the C-Dside plane.
As shown in Figs.10( consistent with the previous stated observations.Overall, the TKE data as captured by LES are considered to be in good agreement with the PIV experimental data, with the mention that the solutions obtained with tripping the boundary layer (LES-WT) approaches better the experimental measurements of the twin square jets.By resolving the finer turbulent flow structures with a higher grid resolution near the wall and jet shear layers, the overestimated TKE level that is observed in the initial jet shear layers can be reduced.This, in turn, can result in changes to the following flow structures, as the LES results can approach the experimental data more closely.In Secs.IV C and IV D, discussions are focused on LES-WT unless there is stated otherwise.The instantaneous acoustic Mach number of the twin jets (u=c 0 , where c 0 is a speed of sound at the ambient air state) and pressure fluctuations in the near-field LES region are presented in Fig. 11.Smallscale turbulence induced by shear-layer instability waves closer to the nozzle and the transition to turbulent state can be observed in the jest shear layers.As convected in the streamwise direction, the developed fluid flow structures grow up to larger scales enhancing the mixing effect between the jets and the ambient air.Power spectral density (PSD) based on the axial velocity component acquired in the lateral jet shear layers at different downstream locations is presented and analyzed in Fig. 12.Note that the frequency is normalized by the ideally expanded jet diameter (D j ) and jet velocity (u j ) (i.e., Strouhal number, St ¼ fD j =u j ).As St increases, the amplitude of PSDs gradually decreases with the slope of around St À5=3 , which follows Kolmogorov's law of turbulent energy cascade.One can find that every spectrum demonstrates rapid drop over a certain high St.It indicates a limitation of the resolution caused by the corresponding grid size and the maximum resolving performance of numerical schemes.
Nevertheless, all spectra follow the À5/3 slope for about one order of magnitude (frequency range), indicating that the grid resolution with the current LES is sufficient for resolving a large range of scales (the most energetic ones) in the inertial subrange.The peak can be found at the St of around 0.37 in the closest location of the nozzle.It is probably due to coherent structures of the axial velocity in the jet shear layer at the screech frequency, which will be continuously discussed in Secs.IV C and IV D.
Returning to Fig. 11, one can find upstream-propagating low frequency waves with a constant wavelength shown in the dashed-box A, representing a generation of screech noise.On the other hand, complicated waves of broadband shock-associated noise (BBSAN) can be observed in the side-line direction.In the flat-side plane (more precisely in the dashed-box B), however, a particular pattern appears.In the downstream regions (for z/h > 10), strong acoustic waves characterized by larger wave lengths (of lower frequency) associated with the large-scale turbulent mixing noise component are visible.
To estimate the development of pressure fluctuations associated with jet-to-jet interactions and the influence of the turbulence tripping on this, a time variation of the pressure fluctuations at the center of the twin nozzle (x/h ¼ y/h ¼ z/h ¼ 0) is investigated for both LES-NT and LES-WT (Fig. 13).It should be noted that the zero in the time axis represents the starting time of data acquisition when it is judged that the jets are fully developed (i.e., 300h/u j ).It can be observed that the amplitudes of the pressure fluctuations for LES-NT significantly vary in time whereas LES-WT case represents more uniform amplitudes.
To examine whether such dominant variations are relevant to a particular frequency or not, a wavelet transform is applied to the pressure fluctuations signals as presented in Fig. 14.The dimensionless frequency of about St of 0.37 is captured, which is assumed to be screech.A distinct mode in the fluid flow can be expected to appear for LES-WT based on the intensive and consistent signals in the wavelet transform analysis.
The contribution of the pressure fluctuations in the near-field can be examined in overall sound pressure level (OASPL), as presented in Fig. 15.It is noted that root mean square of pressure at every resolved frequency was used to generate the OASPL fields, and hydrodynamic as well as acoustic pressure fluctuations are included in the contour levels.According to the contour levels in the flat-and C-Dside planes, the significant larger OASPL levels of around 5 dB are found in the upstream region corresponding to the flat-side plane, region characterized by stronger radiation of upstream propagating waves (i.e., screech noise).In addition, a clear directivity of pressure fluctuations can be found in the side-line direction of twin jets, in agreement with the observed acoustic waves in the dashed-box B of Fig. 11(a).Thus, although further investigation is needed, it can be presumed that the screeching phenomenon is mainly associated with the physics presented in the flat-side plane.It is also interesting that a standing wave pattern is observed between the twin jets (z/h > 0) and in the cavity zone (z/h < 0) as well, which can be generated by superposition of oppositely moving waves.In the present case, this pattern between the twin jets is formed by the interaction of downstreamtraveling instability waves in the jet shear layer (jet shear-layer instabilities) and the upstream-propagating acoustic waves in the subsonic region between the twin jets.On the other hand, the standing waves observed in the cavity zone are estimated to be generated by upstream-propagating acoustic waves and the ones reflected from the nozzle surface located at around z/h ¼ À2.75.Since the reflected waves have the similar amplitude with the original waves, the standing wave pattern in the cavity zone appears to be very clear.

Two-point spacetime cross-correlation analysis on the jet shear layers
To investigate upstream-and downstream-propagating components of pressure fluctuations, a two-point space-time crosscorrelation analysis is performed on time-series data that have been acquired at 40 probe points equally distributed in space within 0 < z/h < 15 in each upper, lateral, and inner jet shear layer.The nomenclature used for the developed jet shear layers is given in Fig. 1.The monitoring probe located at z/h ¼ 4.8 is used as the fixed reference point that stands at approximately fourth shock cell, and all the probes positioned along each jet shear layer are targeted for the analysis.Each correlation result obtained from the two-point space-time crosscorrelation analysis at the corresponding axial location of the targeted point is presented in Fig. 16.
First of all, the highest correlation of 1.0 can be seen at z/h ¼ 4.8 from all the figures because the cross-correlation is locally converted to the auto-correlation at this position.The positive slopes captured in the lateral and inner jet shear layers, of which value is almost the speed of sound (i.e., z/s $340 m/s) at the ambient air state, suggest the presence of upstream-propagating pressure waves in the flat-side plane.They are relevant to a tonal noise with a periodic pattern appearing along the time-lag axis.It seems that the correlation levels of the positive slopes obtained in the inner jet shear layer are higher than those obtained in the lateral jet shear layer.The cause can be associated with the cumulative effects of upstream propagating waves from both the twin jets.On the other hand, negative slopes are also found approximately at downstream locations greater than z/h > 6, which indicates the presence of downstream-traveling pressure waves.Since the periodicity of the negative slopes is the same with the positive ones, the downstream-traveling pressure waves are relevant to the convection of the coherent flow structures in the jet shear layers representing the same frequency of the screech.
On the other hand, it is also important to determine whether the upstream-propagating waves generated by the twin jets are phaselocked or not.Therefore, the two-point space-time cross-correlation analysis is performed on two targeted points that are aligned at the same axial location but placed in each upper shear layers of Jet1 and Jet2.The same can be done for the inner shear layers as well.According to the results representing the unchanged signs of strong correlations along the z/h axis in the inner jet shear layers, it can be estimated that the upstream-propagating waves for LES-WT are strongly phased-locked in the flat-side plane.In addition, it can be also expected that the distribution of pressure waves is symmetric with respect to the symmetry plane (x/h ¼ 0) because the positive correlation (in red) is observed at the time-lag equal to zero.One cycle of the periodic pattern captured along the time-lag axis in Fig. 17 Consequently, based on the cross-correlation analyses, it is expected that the present twin square jets represent a symmetrical flapping mode in the flat-side plane.This explains the symmetrical distribution of upstream-propagating waves and the significantly larger amplitudes of the standing waves shown in Fig. 11(a) and Figs.15(a) and 15(b).The present result is interesting because the opposite behaviors of twin jets have been observed for the twin rectangular nozzle configuration (AR of 2) under the same jet operating condition 42 (i.e., anti-symmetrical flapping mode observed in the C-D-side plane).In addition, the twin square jets experience a comparatively continuous coupling process of the twin jets, whereas the same process in the case of the twin rectangular jets was intermittent. 42Due to the symmetric nature of the instability propagation, the shear layers are coupled in a self-energizing excitation loop leading to continuous production of screech tones.The nozzle geometry is believed to be a significant factor in determining the shear layer along which the instability propagates although the exact mechanism needs further investigation.

Spectral proper orthogonal decomposition (SPOD) analysis
A symmetrical flapping mode with respect to the symmetry plane has been observed for the twin jets.It was analyzed by the signal databased analysis.To reveal and visualize the dominant modes of the twin jets, a spectral proper orthogonal decomposition (SPOD) 58 is additionally performed, which is known as the optimally averaged dynamic mode decomposition (DMD) for providing the best description of the statistical variability of the flow. 58A total of 500 consecutive datasets sampled at equal time-steps of 7:5 Â 10 À6 s are considered for SPOD.The analysis is carried out using first axial velocity data and then pressure fluctuations.
The leading modes of SPOD analysis based on the axial velocity fluctuations are shown in Fig. 18 for the fundamental mode (St $0.37) and for its higher second and third harmonic frequencies.It can be observed that large flow structures dominate the twin jets in the flatside plane at the fundamental frequency.The modes are antisymmetric with respect to each jet centerline, but are symmetric with respect to the symmetry plane, which confirms the suggestions discussed in Sec.IV C 2. In the C-D-side plane, the flow structures seem to be symmetric with respect to the jet centerline.However, one may consider that parts of the large flow structures caused by flapping can appear in the C-D-side plane.The dominant flow structures of the second and third harmonic components are more concentrated in the jet shear layers that include finer scales.Figure 19 presents the SPOD analysis based on the pressure fluctuations for the fundamental frequency and the higher harmonics.The flapping mode is captured; however, more evident structures appear in the leading mode of SPOD of the pressure fluctuations as compared to the corresponding analysis carried out using axial velocity data.It suggests that the pressure field is much subordinated to the screeching phenomenon.Acoustic wave propagations of the second and third harmonic components also can be observed.It seems that the harmonic components have particular directivities in the flat-side plane [see Fig. 19(c)]: the second harmonic wave propagates at an angle of about 30 with respect to the upstream direction, and the third harmonic wave propagates over wider area of upstream and side-line directions.Here, it is again confirmed the responsibility of the particular directivity observed in the OASPL contours originating from the downstream location of 6 < z/h < 8 in Fig. 15(a) is due to the second harmonic waves.It can be also found that the harmonics generate symmetric mode distribution with respect to the symmetric plane.
It was discussed that there are upstream-and downstreamtraveling pressure waves that form the standing waves in the jet shear layers and surrounding area, and both are inevitably included in the leading mode of SPOD.Furthermore, we separate them into upstream and downstream propagating components in the flat-side plane by applying a Fourier decomposition in space.In Fig. 20(a) representing only the upstream-propagating leading modes of SPOD of the pressure fluctuations at the fundamental frequency, it can be found that the large amplitude of waves exist both inside and outside of the jets.The feedback-loop mechanism of screech tone has been accepted as being relevant to an external mode; however, the present twin square jets present a strong internal mode (confined inside the jets) with respect to each jet centerline, as also discussed in recent LES study of the high-temperature single rectangular jet. 32On the other hand, the downstream-traveling waves that include both the coherent flow structures and the reflected acoustic waves from the nozzle cavity, shown in Fig. 20(b), seem to be weaker.
Figure 21 presents the amplitudes of the upstream-propagating and downstream-traveling leading modes of SPOD at the fundamental frequency, which are extracted along the line between jet center and  lateral jet shear layer.In addition, the bounds of shock-cells are also indicated by dashed lines corresponding to the valleys of the timeaveraged axial velocity (see Fig. 7).It can be found that the upstreampropagating wave initiates at approximately z/h ¼ 7.5 that corresponds to the ends of the positive slopes of the cross-correlations shown in Figs.16(b) and 16(c).However, it grows as propagates upstream, and reaches the maximum amplitude near the fifth shock-cell believed to have the largest contribution to the screech phenomenon.It is estimated that the wavelengths of the upstream-propagating mode are strongly relevant to the length of the shock-cell structures, whereas the downstream-propagating mode induced by the instability waves and the coherent flow structures represents a regular wavelength.

D. Noise characteristics of twin square jets
In the experimental facility of University of Cincinnati, far-field acoustics were measured by the microphones located at the radial distance of 67.2h (44 in.) from the center of the twin nozzle, and the acquired data as function of directivity were utilized for the validation of the LES data and Ffowcs Williams and Hawkings (FW-H) calculations.The detail of the FW-H (permeable) surface can be shown in Fig. 11.The size and location of the surface were determined by referring to the previous LES study on single rectangular jets. 29The finest grid size on the FW-H surface is approximately 0.024h, corresponding to St of around 3.5 based on the resolving performance of the present solver (approximately 11 grid cells per wavelength required for propagating acoustic wave 29 ).However, the maximum resolving frequency of far-field spectra will be lower due to varying grid size on the FW-H.To obtain converged (stable) noise spectra at low frequencies, approaches accounting for several closed disk surfaces 59 or corrections of the quadrupole term in the FW-H equation 60 can be considered.However, the classical approach that uses an open disk to reduce computational costs is applied here. 29Time-dependent flow and acoustic data were saved every 100 time step (i.e., 5 Â10 À6 sÞ at all the control points of the FW-H surface, and 8000 consecutive time-series data were utilized for the far-field noise calculation. Figure 22 compares the calculated OASPLs with experiments in far-field (at the radial distance of 67.2h from the center of the twin nozzle) as a function of directivity angles (h < 90 are the upstream angles with respect to the twin nozzle exit plane, while the downstream aft angles are indicated by h > 90).It can be observed that the LES/ FW-H data are in good agreement with the measurements.In both flat-side and symmetry planes, the maximum noise level can be found at around 140 where the turbulent mixing noise can be measured.In particular, one can find the local hump of OASPL variation near 90 in the flat-side plane in both LES/FW-H and the experimental results.
The noise spectra at far-field for the computation and the experiment are presented in Fig. 23.To generate the spectra, Fourier transformations were applied to sub-set data of 250 000 time steps filtered by the Hanning window function.75% overlapping was applied for the sub-set results when the spectra were averaged.On the other hand, acoustic pressure data for 1-s period with a sampling rate of 204 800 were used for generating the spectra in the experimental framework.The same frequency resolution has been targeted (i.e., the same time-range of the sub-set data) in both data sets.The results show that the solutions of the LES/FW-H calculation are accurate in terms of the variations of broadband noise spectra shown at different directivity angles.The overall shapes of the broadband noise spectra obtained in the flat-side plane and the symmetry plane seem to be similar, except the distinct more pronounced hill (shown within 0.5 < St < 1.0) of BBSAN at 90 in the flat-side plane, which may explain the higher OASPL at this directivity angle in Fig. 22(a).Regarding the tonalities, significantly intense tones are observed at the fundamental and the second harmonic frequencies.The spectrum at 90 obtained by LES/FW-H only registers the intense second harmonic tone in the flat-side plane, which could be also observed in the LES result of the twin rectangular jets. 42Otherwise, the major characteristics of the noise are well predicted by LES/FW-H, showing the contributions of the tonal and broadband noise generated by the interaction of twin square jets.
In the flat-side plane where we estimated that the symmetrical coupling mode occurs, the levels of the fundamental and the second harmonic tones are higher than the same components observed in the other plane.In addition, a third harmonic also seems to be captured at upstream angles.The details of the screech tone frequency are given in Table II, which also contains the one predicted by Tam's model 16 suggested for a rectangular jet as where b and h are width and height of the nozzle exit, respectively (b ¼ h for a square nozzle exit).a 1 is the speed of sound at the ambient state of air, and u c is convection velocity of turbulence in the jet shear layers, which is assumed as 0.79 u j based on our previous study. 11It can be found that the screech frequencies show good agreement even though the result obtained by the experiment presents a slightly higher value.
In Fig. 24, the far-field noise spectra obtained at continuous directivity angles are shown to present the noise characteristics of the twin square jets.In both the flat-side and symmetry planes, prominent contributions of the large-scale turbulent mixing noise-the large amplitude at very low St and at large downstream/aft directivity angles above 130 -can be found.In addition, the BBSAN-the large amplitude initiated around St of 0.5 at the low directivity angles and gradually shifted toward higher St as the directivity angle increases-is also observable.In addition, the fundamental screeching tone and the second harmonic can be clearly visible along the directivity angle axis in both reference planes.However, they are stronger in the flat-side plane with introducing the third harmonic as well.The tonalities seem to have different amplitude levels and to be partially visible.Thus, the fundamental tone and the second harmonic component are separately plotted as shown in Fig. 25.Note that a moving average filter was adopted to make plot smoother.The results show that the fundamental tone observed in the flat-side plane has low amplitude near 90 degrees side angle directivity and has a larger amplitude when moving at lower (upstream with respect to the nozzle exit plane) or higher (downstream of the nozzle exit plane) directivity angles.In the symmetry plane, the amplitude of the fundamental tone is rapidly reduced upstream, whereas it increases at the wider angles downstream (consistent with the SPOD of pressure fluctuations presented in Fig. 19).On the other hand, the opposite trend of the directivity can be found in the results of the second harmonic component.

V. CONCLUSIONS
In this paper, the aerodynamics and aeroacoustics of the twin square jets were investigated under the screeching jet conditions (NPR of 3.0 and TR of 1.0).To better represent the disturbed jet shear layers induced by the turbulent boundary layers at the walls inside the nozzle, which is believed to be the situation in the experiment, a small step was applied in the convergent section of the nozzle in the wall boundary layer for generating flow disturbances.As compared to LES data without tripping (LES-NT), by tripping the boundary layer inside of the nozzle (LES-WT), the fluid flow solution approaches better to the experimental data in terms of the length of the supersonic region that is influenced by the shock-cell structures and TKE levels in the jet shear layers.Qualitative and quantitative analysis of the captured screech noise phenomenon has been carried out by analyzing both snapshots (the Mach number and the pressure fluctuations fields) and by performing statistical analyses (frequency spectra and mode decomposition).It has been found that the screech noise is more active in the flat-side plane.The contribution of the pressure fluctuations in the near-field depending on the observation planes was examined by OASPL contours.In addition, the high amplitudes of the standing waves were observed in the region between the twin jets and in the cavity region shaped by the outer surface of the nozzle.The results of the two-point space-time cross-correlation analysis showed that the screeching phenomenon mainly occurs in the flat-side plane, and the coupling mode of the twin jets was determined as the symmetrical flapping mode.Through the SPOD analysis, we confirmed the symmetrical flapping mode of the twin square jets confined inside jets and visualized the flow structures and the directivities of near-field acoustics at the fundamental and higher harmonic frequencies.The downstream-traveling modes depicted a regular wavelength along the stream-wise direction, while the upstream-propagating mode is connected to the length of the shock-cell structures.Satisfactory validation was achieved through comparison with the measurement data.The screech tone has higher amplitudes in the flat-side plane where the main coupling of the twin jets was determined.Around 90 , the fundamental screeching tone becomes weak and is increased at the lower  or higher directivity angles, whereas the second harmonic component presented an opposite trend.

FIG. 4 .
FIG. 3. Nozzle details in the C-D-side plane (y-z plane) showing the step used for tripping the boundary layer: (a) the location of the step inside of the nozzle; and (b) step geometry and computational grid details around the step.

FIG. 5 .
FIG. 5. Turbulent intensity at the nozzle exit plane in the (a) flat-side plane and (b) C-D-side plane.The blue and red lines indicate LES-NT and LES-WT cases, respectively.

FIG. 8 .
FIG. 8. Comparisons of time-averaged axial velocity profiles at several axial locations in the (a) flat-side and (b) C-D-side planes; LES-NT vs LES-WT vs PIV.

FIG. 10 .
FIG. 10.Comparisons of normalized turbulent kinetic energy (TKE) profiles at several axial locations in the (a) flat-side and (b) C-D-side planes; LES-NT vs LES-WT vs PIV.

FIG. 11 .
FIG. 11.Instantaneous acoustic Mach number and over-imposed pressure fluctuation fields obtained by LES-WT in the (a) flat-side and (b) C-D-side planes.
(b) is the same with the previous results of Figs. 15 and 16.

FIG. 14 .
FIG. 14. Wavelet transformations of the pressure fluctuations at the center of twin nozzle (x/h ¼ y/h ¼ z/h ¼ 0) for (a) LES-NT and (b) LES-WT.

FIG. 16 .
FIG. 16.The cross-correlation analysis on pressure fluctuations along the (a) upper, (b) lateral, and (c) inner jet shear layers for LES-WT [see also Fig. 1(b) for the nomenclature used to define the jet shear layers].

FIG. 17 .
FIG. 17.The cross-correlation analysis on pressure fluctuations along the (a) upper jet shear layers and (b) inner jet shear layers of the two jets for LES-WT.

FIG. 18 .
FIG. 18.The leading mode of SPOD of the axial velocity fluctuations at the (a) and (b) fundamental, (c) and (d) second harmonic, and (e) and (f) third harmonic frequency in the (a), (c), and (e) flat-side and (b), (d), and (f) C-D-side planes.

FIG. 19 .
FIG. 19.The leading mode of SPOD of the pressure fluctuations at the (a) and (b) fundamental, (c) and (d) second harmonic, and (e) and (f) third harmonic frequency in the (a), (c), and (e) flat-side and (b), (d), and (f) C-D-side planes.

FIG. 21 .
FIG. 21.The amplitudes of the upstream-propagating and downstream-traveling leading modes of SPOD of the pressure fluctuations at the fundamental frequency.The plots are extracted along the line between jet center and lateral jet shear layer (i.e., x/h ¼ 1.31).

FIG. 23 .
FIG. 23.Noise spectra at far-field (67.2h) for different directivity angles; LES/FW-H as compared with the experimental microphone data in the (a) flat-side plane and (b) symmetry plane.The scale exacts at 45 and a 45 dB offset is applied between each angle.

TABLE I .
Jet operating conditions.M j ; u j ; and Re indicates the ideally expanded Mach number, the fully expanded jet velocity, and the Reynolds number based on the equivalent diameter (D eq Þ, respectively.

TABLE II .
Screech tone frequency; comparisons between the Tam's model, experimental data, and LES-WT data.