Unsteady shock front waviness in shock-buffet of transonic aircraft

The unsteady transonic aerodynamics of a wing-body configuration is investigated by solving the Unsteady Reynolds-averaged Navier-Stokes equations closed with the full Reynolds Stress Model. This work presents the prediction of flow field characteristics during deep shock-buffet penetration of a transport aircraft-representative geometry. Mach number of 0.85, and Reynolds number of 5 million based on the mean chord, are selected to reproduce experimental test conditions that serve as validating datasets. The results obtained give information about both surface and flow field shock-buffet dynamics. An unsteady shock front is observed on the suction side of the wing which gives birth to the so-called buffet cells. Flow field characteristics are dominated by the presence of lambda-shaped shocks and fully separated boundary layer over a significant part of the wing.

With advancing technology, the use of Computational Fluid Dynamics (CFD) in shock-buffet prediction has increased substantially in the last two decades. The 2D experimental results for the OAT15A aerofoil [8] have been used for validation in a series of CFD studies in which the Unsteady Reynolds-averagedNavier-Stokes (URANS), Large Eddy Simulations (LES), and hybrid RANS-LES approaches have been applied [9][10][11][12][13][14][15][16]. The variation in results presented revealed that there is still an uncertainty in prediction of shock-buffet by means of URANS, particularly when closed with eddyviscosity-based turbulence models. CFD has also been successfully applied to shockbuffet prediction for 3D configurations. Iovnovich and Raveh [17], produced URANS results which show the effect of varying sweep and aspect ratio of rectangular wings on shock-buffet. They also coined the term of buffet cell which is now used to describe the spanwise extent and oscillation of the shock wave. The presence of these spanwise oscillations have since been observed in detail by Sugioka et. al. [18] in a transport aircraft-representative configuration.
A wing-body aircraft configuration was simulated in shock-buffet by Sartor and Timme [19,20] using both URANS and Delayed Detached Eddy Simulations (DDES). Both approaches compared well with experiments. Recently, the NASA Common Research Model (CRM) 1 , representative of modern wide-body aircraft, was used as a test bed and simulated using a Zonal DES approach [21][22][23]. The three studies which compare numerical results with experimental observations of Sugioka et. al. [18] and Koike et. al. [24] showed intricate shock-buffet dynamics and a qualitative comparison with experiments. Similar observations were presented by Ribeiro et. al. [25] who used the Lattice Boltzman Method (LBM). Timme [26] produced a shock-buffet onset prediction on the NASA CRM which was in good agreement with experiments. He also found evidences that reinforce Crouch's global instability shock-buffet onset mechanism. Still no efforts have been made to further characterise this phenomenon evolution with increase in flight angle of attack.
The current contribution introduces URANS-based simulations in the context of predicting shock-buffet dynamics on industrially-relevant aircraft configurations. A second-moment closure, the full Reynolds Stress Model (RSM), is selected as a turbulent closure for the URANS equations. RSM are considered to be part of an advanced class of turbulence models, not widely used due to their cost. However, in this case, the use of RSM is justified. Gianellis et. al. [15] found that RSM closures are more robust than typically used Spalart-Allmaras (SA) and Shear Stress Transport (SST) turbulence models for 2D shock-buffet prediction. Steady-state simulations by Apetrei et. al. [27] and many participants at the American Institute of Aeronautics and Astronautics CFD Drag Prediction Workshop 2 (AIAA CFD DPW) found that for this particular test case, off-design flight predictions produced by closing the RANS equations with eddy-viscosity turbulence models can result in unphysical side-of-body separation [28].
Experimental campaigns of National Aeronautics and Space Administration (NASA) [29,30] provide the base for validation of current results, including the decision on flight conditions. The observations from Japan Aerospace Exploration Agency (JAXA) wind tunnel campaigns [18,24] are used for qualitative comparison, since the difference in flight Re and wind tunnel far field dynamic pressure resulted in different wing deformations, making a quantitative comparison not possible [31].
The remainder of the paper is structured as follows. Section 2 introduces the details of the numerical approach followed by a presentation and discussion of results in Section 3. Conclusions are then drawn in Section 4.

Geometry
The wing-body NASA CRM [32] from the 6th AIAA CFD DPW, shown in Fig. 1, is used as a representative configuration of modern transport aircraft. For the current study, the geometry with applied wing deformations measured in the European Transonic Wind tunnel (ETW) at 4°incidence is used to improve the quality of the comparison between numerical and wind tunnel data [33]. The full dimensions of the aircraft, given in Table 1, reveal its resemblance with today's wide-body jets.
Experimental data documenting the shock-buffet regime on this aircraft is available from Balakrishna and Acheson [30] for M and mean-chord Re of 0.85 and 5 million respectively. A smaller wind tunnel model was also experimentally tested in JAXA wind tunnels at lower mean-chord Re of 0.96 and 1.5 million. Although the latter wind tunnel campaign is more thorough, only a qualitative comparison can be produced since wind tunnel model deformations differ [31].

Grid
The grid used in this study has 17.5 million nodes and is produced using the commercial meshing tool ANSYS ICEM-CFD. The half-modelwing-body aircraft is placed at the centre of a semi-spherical fluid domain with a radius of 100 semispans. A multiblock structured grid is generated using quad and hex elements only in order to maintain mesh quality for viscous flow. The average y + value near wall is maintained to be lower than 1 as required by the implementation of the mathematical approach used [34]. This grid is a modified version from a previous study published by the authors which also includes a mesh-sensitivity analysis in line with AIAA CFD DPW guidelines [27]. The surface mesh and its features can be observed in Fig. 2.

Solver
The finite volume unstructured RANS solver DLR-Tau [35] is used in conjunction with the DLR-FlowSimulator [36] python interface. The convective fluxes are discretized using a second-order upwind scheme. Time-marching simulations are run with a LU-SGSsemi-implicit time stepping scheme.
The RANS equations are closed with the SSG/LRR-g variation of the RSM model developed at DLR [34]. The model is a blend of the Speziale-Sarkar-Gatski (SSG) [37] and Launder-Reece-Rodi (LRR) [38] RSM turbulence models, with LRR-like behavior near-wall, and SSG-like behavior away from the wall [39]. The blending of the two models, originally in the SSG/LRR-ω model, has shown promising prediction of transonic aerodynamics in the DLR contributions to the AIAA CFD DPW [40,41]. The current model solves one equation for each Reynolds stress component, and an additional equation for a turbulent scale. In this case, the turbulent scale is g ¼ where ω is the specific rate of dissipation. The transformation to g solves numerical instabilities that the SSG/LRR-ω model had exhibited. This transformation did not have an impact on the accuracy of the results as shown in the work by Togiti and Eisfeld [34]. Initial simulations were also performed by closing the RANS equations with the SA-negative [42] turbulence model.
For the URANS computations, time-marching simulations were performed with a dual time-stepping approach with 100 inner iterations per physical time step. Solution convergence was monitored through the lift coefficient residual, which decreased by at least 6 orders of magnitude. A time step Δt of 0.5 ms was chosen based on the physical convection times of an air molecule travelling over the mean chord at free stream velocity, and on the buffet frequencies expected. The total simulation time for all incidences was at least 1 s. This resulted in at least 3 periods of the lowest frequency identified, and more than 15 of the highest. The time step sensitivity was evaluated by reducing Δt to 0.1 ms and comparing the unsteady buffet response.

Computational test matrix
For the present studies, the Mach number is maintained constant at 0.85. Mean-chord Re of 5 million is set according to the NASA test. The incidence α was increased from 4°to 5°, 6 and finally 6.5°. Steady-state results were obtained at 4°incidence and will be presented in Section 3.1 as a means of validation. Unsteady time-marching simulations were restarted from semi-convergedsteady-state RANS results. Simulations at lower  incidences are used as initialization for higher incidence ones (e.g. the results obtained at α = 4°are used to initialize the run for α = 5°, and so on).

Steady RANS
Steady-state results are produced for the incidence α = 4°, a comparison of the SA-negative [42] and SSG/LRR-g RSM results with experimental observations is possible. Figure 3 presents pressure coefficient (C P ) distributions extracted at 6 spanwise (η) locations. Experimental results were reported by Rivers and Dittberner [29] and were obtained in the NASA wind tunnel facilities; they are easily accessible on the NASA CRM webpage. In a first observation, we re-iterate that the SA-negative produces large discrepancies at the wing root (η = 0.131 and 0.283) due to an overpredicted side-of-body separation. This problem is eliminated by applying a Quadratic Constitutive Relationship (QCR) [43] which transforms the SA model into a non-linear one that is better able to capture the flow physics at the intersection between the wing and the fuselage. The benefits of using QCR for corner flows have been analysed more in depth by Dandois [44], to which the reader is referred.
Even after enabling QCR, the SSG/LRR-g-based results compare better with experiments than those obtained with the SA model. The two inboard slices (η = 0.131 and 0.283) are predicted to nearly identical as experiments. The computed results show stronger shocks being computed in the mid-wing region (η = 0.502 and 0.727), however the pressure recoveries downstream are well predicted. Towards the wing tip, CFD predicts higher loading, visible through a lower pressure plateau. At this incidence, Balakrishna and Acheson [30] reported small amplitude pressure fluctuations, however, the current simulations did not capture any unsteadiness after switching to URANS mode. It is currently believed that the discrepancies can be attributed to unsteadiness due to shock-buffet onset which is not modelled by our current mathematical and numerical approaches.

Time step-sensitivity analysis
Time step sensitivity is assessed via time-marching simulations at M = 0.85, Re = 5.0 million, and α = 6.0°. Simulated flow field data output at these conditions is used as  Figure 4 shows the extracted results at the spanwise locations η = 0.5 and 0.6. On the left hand side, each line corresponds to the time-accurate C P recorded at a chordwise location. On the right hand side, the root mean square (RMS) of fluctuation, C Prms ' = rms (C P,mean -C P ), of those signals is given in comparison. The amplitude of oscillation is observed to vary locally but does not change the overall phenomenon. The location and frequency of maximum oscillation remain consistent. In a trade-off between computational power required and benefits of a smaller time step, it was decided that for this study, Δt = 0.5 ms was sufficient for the purpose of this study.

Global forces and surface sampling
Balakrishna and Acheson [30], measured buffet-onset on the CRM through small amplitude pressure/strain fluctuations at an incidence of 3°. A deep-buffet penetration regime, characterized by large amplitude pressure fluctuation followed as the angle of attack passed 5°. The characteristics of shock-buffet onset were previously numerically predicted by Timme [26]. This study focuses on the second regime, where shock-buffet is expected to be fully developed (α > 5°).
The computed lift coefficient, C L and comparison with the mean C L from experiments are given in Fig. 5 for incidences α = 5, 6, and 6.5°. Figure 5a shows the evolution of the time-accurate C L at all three incidences. The first 0.2-0.3 s of simulated time reveal a transition between the partly-converged RANS and URANS solutions. It indicates that the RANS solution is not the mean of the URANS one. At all three incidences, the RANS solution fluctuated around a lower C L than URANS. Secondly, the unsteadiness captured in Fig. 5a is not harmonic as previously seen for the RBC12 aircraft [19] or in 2D studies for the OAT15A aerofoil [9][10][11][12][13][14][15][16]. Although not shown here, the drag coefficient, C D , and pitching moment coefficient, C M , exhibit similar patterns. There is a clear unsteadiness in the global forces, however this is not a conclusive evidence of shock-buffet. On the right hand side, Fig. 5b gives the mean C L and RMS of fluctuation (error bars) in comparison with mean experimental values from NASA wind tunnel campaigns. The computed C L compares well at an incidence of 5°, however at higher incidence the values are over predicted by 1-2 lift counts (1 lift count = 0.01 in C L ). The error bars, indicating C LRMS ' show that with an increase in incidence, the amplitude of fluctuation in C L increases too.
A better detection of shock-buffet is achieved by sampling surface pressure to identify the areas where flow unsteadiness is present. On-the-fly pressure probing is performed. The DLR-FlowSimulator python interface allows for the definition of 64 surface probes (see Fig. 6) placed at the location of the unsteady pressure transducers in the wind tunnel experiments [18,24,30]. C P is extracted every 5 timesteps resulting in an extraction frequency of 400 Hz, equivalent to a non-dimensional Strouhal number, St = fc mac /U ∞ , of 9.5; sufficient to collect 20 samples per period of oscillation at highest expected frequency based on previous experimental observations. The surface probes are inside the region of η = 0.5-0.731 and local x/c = 0.16-0.8, and benefit the current work by giving To conclude the surface analysis for shock-buffet detection, computed C P slices are presented at 6 spanwise locations against mean experimental data in Fig. 10. It is visible that at η = 0.727-0.846 sections the shock location is predicted marginally more upstream than in the wind tunnel tests. The pressure sides, suction peaks and pressure recoveries downstream of the shock are overall well predicted. Starting with η = 0.397, the timeaccurate C P reveals significant unsteadiness in the region downstream of the shock. In the mid-wing region of η = 0.502-0.603, this unsteadiness is coupled with significant shock travel and a variation in trailing-edge loading on the pressure side of the wing. Outboard at, η = 0.727-0.95, the unsteadiness reduces to mainly an oscillation of the shock.

Shock-buffet dynamics
Because shock-buffet data at current flight conditions and with current CRM configuration is limited, this section focuses on analyzing shock-buffet characteristics by concentrating on shock dynamics on the wing surface and inside the flow field. The results computed for α = 6°are presented herein. This paper comes accompanied with supplementary animations that help visualise flow unsteadiness. They are available on the Cambridge Journals and Sheffield Online Research Data (ORDA) repository webpages, DOI link available [45].   [18,24] and NASA [30] plus additional ones included in CFD Time-accurate C P is not reported herein, however an animation is provided as supplementary material. It shows that the shock oscillates on the wing surface in a sinusoidal manner. To better visualize this oscillation, buffet cells contours are provided in Fig. 12. Contours of C P ' = C P -C Pmean of alternate colours indicate if the local shockposition is upstream (+ve C P ') or downstream (−ve C P ') when compared to the mean location. The computed buffet cells are qualitatively similar with the wind tunnel observations of Sugioka et. al. [18] who have experimentally tested the CRM at a lower Re. Figure 12 shows that each buffet cell travels towards the wingtip; a black * marks the most inboard buffet cell and aids the reader observe its travel in time. The distance between two buffet cells of same colour is then the wavelength of oscillation in the spanwise direction, λ, measured as 0.2 semispan. Knowing the frequency of oscillation which was identified previously in Fig. 8b St = 0.5-0.6 then the convection speed of the shock, U C = λf can be approximated as 0.4 to 0.5 of the freestream velocity, U ∞ .
The flow field output is analysed to determine the shock-buffet dynamics in the flow field. Figure 13 shows an instantaneous snapshot of the shock, extracted using the normal Mach number criteria developed by Lovely and Haimes [46]. The top view in Fig. 13b shows that the shock extends over the entire span of the wing. Inboard of the trailing-edge break, the top view of the shock-front shows that the main shock is preceded by a weaker one. From the top the shock front appears to be wide however this is because in the mid-wing region, the shock structure is λ-shaped near the wing surface, becoming normal only far away. This is visible in Fig. 14 where slices of the shock isosurface are plotted at η between 0.283 and 0.95. We appreciate that the shock extraction algorithm is not perfect, however the lambda shock structure near the wall is easily visible at η = 0.283 and 0.397. Close to the wingtip, at η = 0.95 the shock is normal and remains as such in time.
An animated version of Fig. 14 is provided with the supplementary material. Figure 15 presents two snapshots superimposed. The shock unsteadiness in the flow field is in good agreement with what was observed on the surface via C PRMS ' and C P ' contours in Figs. 11 and 12. Insignificant shock displacement is present at η = 0.283, the region where it is most prominent is η = 0.502-0.727. The out-ofphase oscillation is also visible, in particular if the reader refers to Fig. 15b, c, d and e where local shocks are observed to oscillate in anti-phase. At t = 1.22 s, the shocks at η = 0.397 and 0.603 are downstream, whereas at η = 0.502 and 0.727 are upstream when compared to their locations at t = 1.2 s.
We now address the topic of shock-induced separation and its characteristics. In the 2D studies based on the OAT15A aerofoil, the boundary layer behind the shock presented a periodic separation & reattachment behaviour. However, for this configuration, the boundary layer remains fully separated with complex periodic vortex shedding spanning a significant portion of the wing. Figure 16 gives instantaneous isosurfaces of Q-criterion (the second invariant of the velocity gradient tensor) coloured by Mach number. An animated version of this is also included as supplementary material. The wingtip vortex and additionally another vortex at the trailing edge break are present in addition to successive vortical structures present in between.

Conclusions
The current work assessed the prediction of flow field characteristics during deep penetration into the shock-buffet regime of the NASA CRM aircraft. Mach number and mean chord Re were set as 0.85 and 5 million respectively; incidence was varied between α = 5 and 6.5°. Simulations were produced by solving the URANS equations closed with a second moment closure. Validation of the results showed quantitative agreement with mean experimental data and compared well qualitatively with shockbuffet dynamics observed at lower Reynolds numbers. Time-accurate surface pressure coefficients were extracted and used as on-theflyshock-buffet detection technique. The analysis of the signals revealed that high amplitude fluctuations, with multiple dominant frequencies between St = 0.05 and 0.6, were present at shock location; the frequencies were within the expected range of shockbuffet phenomenon. The surface shock-front was globally unsteady, exhibiting periodic patterns in the spanwise direction that gave birth to buffet cells. Flow field analysis, including the extraction of the shock surface by means of Normal Mach criterion provided a new insight into the flow field dynamics. The shock strength decreased from root to tip; it contained λ-shaped foot over a significant part of the wing. Timeaccurate snapshots showed how the buffet cells are reflected in the flow field. Locally the shock oscillation was out-of-phase between different spanwise locations.
Compared to previously reported results for other geometries, this configuration exhibited fully separated boundary-layer that did not reattach during the buffet cycle. Successive shedding of vortical structures were identified by means of Q-criterion isosurfaces.

Nomenclature
Symbols α = angle of attack,°Δ t = time step, s η = non-dimensional spanwise location μ = free stream viscosity, kg/ms ρ = free stream density, kg/m 3 λ = wavelength, m a = speed of sound, m/s b/2 = semi-span, m c = local chord, m c mac = mean aerodynamic chord, m C D = drag coefficient C L = lift coefficient C M = pitching moment coefficient C P = pressure coefficient f = frequency of oscillation, Hz l = reference length, m Re = ρc mac U ∞ /μ = mean-chord Reynolds Number U C = convection speed U ∞ = freestream velocity, m/s St = f c mac /U ∞ = Strouhal number x/c = non-dimensional chordwise location y + = non-dimensional first cell height Sub−/superscripts RMS = root mean square mean= mean value '= fluctuating component