12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

TURBULENT TRANSONIC BUFFET ONSET PREDICTION ON THE NASA
COMMON RESEARCH MODEL VIA GLOBAL STABILITY ANALYSIS

Andrea Sansica
Aircraft Lifecycle Innovation Hub, Digital Integrated Design Technology Team
Aviation Technology Directorate
Japan Aerospace Exploration Agency
Chofu Aerospace Center, Chofu, Tokyo, 182-8522 Japansansica.andrea@jaxa.jp

Atsushi Hashimoto
Aircraft Lifecycle Innovation Hub, Digital Integrated Design Technology Team
Aviation Technology Directorate
Japan Aerospace Exploration Agency
Chofu Aerospace Center, Chofu, Tokyo, 182-8522 Japanhashimoto.atsushi@jaxa.jp

ABSTRACTerning mechanisms remain unclear. Numerical simulations
particularly suffer of the complexity associated to the buffet
phenomenon and computationally inexpensive tools are constantly sought for its characterization and prediction.
Global stability analysis (GSA) has recently proved agood alternative to classical unsteady RANS (URANS) simulations for buffet prediction (Crouch et al., 2009, 2019; Paladini et al., 2019; Sansica et al., 2022). However, due to thehigh RAM memory costs associated with the resolution of thelinearized dynamical system eigenvalue problem, the configurations studied were simplified to 2D aerofoils or periodic
extruded wings. Timme (2020) recently performed a fully-3D
GSA study on the NASA Common Research Model (CRM),
but only two angles of attack (AoAs) were investigated on an
8 million cells grid. While a 3D buffet-cells mode was foundto be unstable, no evidence of a 2D shock-oscillation mode
was given.
In the present work, a fully-3D GSA study is carried out
with the intent to predict buffet onset for the NASA CRM
body-wing-tail configuration. Due to the combination of atime-stepper and an Arnoldi method, the RAM requirements
are relaxed and a parametric studied on the AoA is performed
to predict the buffet onset on grid resolutions spanning from 3
to 30 million cells. Buffet is characterized both in the vicinity
and far from the onset and the numerical results are comparedwith the experiments carried out at the Japan Aerospace Exploration Agency (JAXA).

A fully-3D global stability analysis (GSA) is performedon the NASA Common Research Model (CRM) body-wingtail geometry. Differently from other methods present in theliterature, the combination of the solver linearization approachand modal decomposition based on time-stepper and Arnoldi
methods allows to relax the RAM requirements and perform
an angle of attack (AoA) parametric study to determine buffet
onset. The RANS solutions, buffet onset and frequencies predicted by GSA are all in good agreement with the experiments
carried out at JAXA’s 2m×2m Transonic Wind Tunnel. While
only a 3D buffet cells mode is found unstable in the vicinityof the onset, additional modes appear at higher AoAs. Besides a 2D low-frequency mode, another high-frequency unstable mode is localized in the separated region at the wing-bodyjunction. This is the first fully-3D GSA study able to predict
and characterize buffet onset over a large range of AoAs.

INTRODUCTION
For commercial aviation, the interaction between shockwaves and boundary-layers represents a critical aspect for bothdesign and certification phases. A specific type of shock-wave
and turbulent boundary-layer interaction is the so-called turbulent transonic buffet. Turbulent transonic buffet is a threedimensional (3D) phenomenon consisting of two distinct instabilities (Lee, 2001; Giannelis et al., 2017). The first one isrelated to the two-dimensional (2D) streamwise shock oscillations, that cause the shock-wave forming on the wing suction side to move periodically back and forth in the streamwise direction. The second one concerns the 3D cross-flow
propagation of the so-called buffet cells. These instabilities
can origin wing vibrations and represent a limiting factor for
the flight envelope. Due to the coexistence of 3D effects, separated turbulent boundary-layers and unsteady shock-waves,buffet represents a very challenging problem and, despite thepast and recent experimental and numerical studies, the gov-

1

PROBLEM FORMULATIONThe compressible 3D RANS equations for a perfectgas can be written using Boussinesq hypothesis in the nondimensional form as:∂q= N (q),
∂t

1

(1)

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
where q = [ρ, ρu, ρE, ρνt ]T is the state vector in the conservative form (with ρ, u, E and νt being fluid density, velocity
vector, total energy, and kinematic turbulent viscosity, respectively) and t is the time. The differential nonlinear RANS operator N can be explicitly expanded as

By choosing the normal mode or wave solution q0 (x,t) =b(x) exp (λt) + c.c., the eigenproblem L q
b = λqb is obtained.
q
The complex eigenvalue can be split in its real and imaginaryparts λ = σ + iω, where σ is the temporal growth rate and ω
the angular frequency. While the angular frequency characterizes the oscillatory behavior, the temporal growth rate indicates whether the equilibrium state bifurcates to another solution. This bifurcation is expressed in a linear framework by theexistence of eigenmodes with a corresponding positive growth
rate.




ρu


ρu ⊗ u + pI −ττ −ττ RN (q) = − ρEu + pu −ττ u −ττ R u + q + qR µ+ρνtρνt u − σM ∇νt + SM

(2)

with

2
1
2u·u
p = (γ − 1)ρE −
Tτ = µ (∇ ⊗ u + ∇ ⊗ u ) − 23 (∇ · u)I τ R = µt (∇ ⊗ u + ∇ ⊗ uT ) − 32 (∇ · u)IµC
q = − Prp ∇Tµt Cp
qR = − Pr∇Tt

While the calculation of the base flow solutions is carriedout by using a classical RANS (nonlinear) solver, the global
stability analysis requires a URANS linearized solver. The
characteristics and numerical strategies used for both solvers
are described below.

(3)

2.1

1 + Ts
T + Ts(4)∗ , with T ∗ the dimensional free-stream
where Ts = 110.4K/Ti,∞i,∞stagnation temperature (the superscript ∗ indicates dimensional quantities). The array of the streamwise, vertical andtransverse directions is indicated by x = [x, y, z]T . Note that all
variables are Reynolds averaged, except for u and E that are
Favre (density-weighted) averaged. The formulations of thecoefficient σM and the turbulent source terms SM depend on
the turbulence model.

1.1

Nonlinear Solver

As RANS nonlinear solver, JAXA’s unstructured-gridflow solver FaSTAR (Hashimoto et al., 2012; Ishida et al.,2017) is used. The governing equations are the compressible Navier-Stokes equations. The cell-center finite volume
method is used for the discretization. The numerical flux
is computed with the SLAU scheme and a dual-time stepping method (Visbal & Gordnier, 2000) is used to perform
an accurate time calculation with an implicit time integration. The LU-SGS scheme (Sharov & Nakahashi, 1998) isused for the pseudo time sub-iterations and the physical time
derivative is approximated by a three-point backward difference method. The Spalart-Allmaras turbulence model (Spalart& Allmaras, 1992) with rotation correction (SA-R) (DaclesMariani et al., 1995) and quadratic constitutive relation 2000
version (Spalart, 2000) is used to close the averaged Reynolds
stresses. The boundary conditions used are: no-slip velocity and adiabatic temperature on the aicraft walls; far-fieldboundary conditions are employed at the lateral boundaries
and the AoA is applied at the inflow of the numerical domain.
All steady solutions are computed by using a large CourantFriedrichs-Lewy (CFL) number equal to 10. The selective frequency damping (SFD) method (Richez et al., 2016; Akervik
et al., 2006) is used to further converge the steady fixed point
solution and filter any possible unsteadiness.

being p the pressure, τ the stress tensor, τ R the Reynolds stresstensor, C p the heat capacity at constant pressure, µ the dynamic viscosity, µt the eddy viscosity, Pr and Prt the classical and turbulent Prantdl numbers, T the temperature, q theheat flux and qR the flux of diffusion of turbulent enthaply.
The Prandtl numbers are considered constant and equal toPr = 0.72 and Prt = 0.90. The dynamic viscosity is assumed
to follow Sutherland’s law asµ = T 3/2

NUMERICAL METHOD

Stability Problem

The GSA problem is based upon the use of the linearized
RANS equations. The first step to obtain this linearized setof equations is to assume that the nonlinear system in Eq. (1)admits an equilibrium solution, qb , defined by N (qb ) = 0
and referred to as fixed point or base flow. In this case, thesteady RANS solution corresponds to the base flow. The
standard small perturbation technique is used to decompose
the instantaneous flow into base flow and small disturbancesq(x,t) = qb (x) + εq0 (x,t), with ε  1. By assuming that the
perturbations are infinitesimal, all nonlinear fluctuating terms
are ignored and the linearized RANS equations can be writtenas

2.2

Linearized Solver

The expression of L q0 is extremely complicated andwould mean a lengthy implementation and modification of theexisting nonlinear version of FaSTAR. To avoid this, a strategy
based on a finite difference method consists of using the nonlinear solver in a black box manner and approximate L q0 via
repeated evaluation of the residual function N (q). A Taylor
series expansion with a first order approximation allows then

L q0 =∂ q0∂t= L q0 ,1[N (qb + εq0 ) − N (qb )],
ε

(6)

(5)
with ε being a small constant. This method and selection of the
ε constant are further discussed in (Tezuka & Suzuki, 2006;
Mettot et al., 2014). A matrix-free method (Edwards et al.,
1994; Bagheri et al., 2009) is then used to solve the eigenb = λqb. Being L the discrete form of L , it is posproblem L q
sible to introduce the exponential propagator M = exp (L∆T )Twhere q0 = ρ 0 , ρ 0 ub + ρb u0 , ρ 0 Eb + ρb E 0 , ρ 0 νt,b + ρb νt0
isthe state vector of conservative perturbation variables andL = ∂ N /∂ q |qb is the Jacobian operator obtained by linearizing the RANS operator N around the base flow qb .

2

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 1. Experimental PSP and RANS pressure coefficient contours on the main wing (left) and pressure distributions at different
spanwise wing sections for AoA = 4.8◦ (right).

that linearly advances the perturbation solution in time asq0 (t n+1 ) = Mq0 (t n ), with t n+1 = t n + ∆T . An Arnoldi algorithm (Arnoldi, 1951; Lehoucq et al., 1997; Barkley et al.,2008) is coupled to the linearized solver (Loiseau et al., 2014;
Guiho et al., 2016; Sansica et al., 2018) to extract the leading
eigenmodes of M. While all conservative quantities are perturbed, the perturbation eddy viscosity is not included in thelinearized snapshot matrix (Krylov base).
The linearized solver, referred to as FaSTAR-GSA, hasbeen validated against the literature for 2D laminar incompressible cylinders, 2D turbulent transonic buffet and applied
to fully-3D unswept and swept wings in the presence of lateralwalls (Sansica et al., 2022). The combination of the solver linearization approach with the modal decomposition techniqueselected allows to: a) leverage all FaSTAR’s high-performancecomputing features for complex geometries, b) be flexible inthe selection of spatial schemes and turbulence models (opening up to sensitivity studies) and c) investigate high-resolution
numerical grids.
Figure 2.

Lift-coefficient based buffet onset prediction.

FLOW CONDITIONS AND NUMERICAL GRID
The numerical flow conditions are selected to match theexperiments carried out at JAXA’s 2m×2m Transonic Wind
Tunnel (JTWT1) on the NASA CRM at a Mach number of 0.85
and a chord Reynolds number of 2.27×106 by Sugioka et al.
(2021). Buffet onset was experimentally determined by using the lift-coefficient method and wing-strain gauge measurements at AoA ≈ 3.5 and 3.7◦ , respectively. Unsteady pressuresensitive paint (PSP) and pressure taps at different wing sections were used to characterize the pressure fluctuations. Buffet cells were observed in a frequency range of 400 − 850 Hz(or, in nondimensional Strouhal number, St≈0.2-0.5).

RESULTSThe geometry selected is the NASA CRM model in thebody-wing-tail configuration. The main wing considers thedeformation measured experimentally at AoA = 2.94◦ . Coarse(≈ 2.8 × 106 cells), medium (≈ 9 × 106 cells) and fine (≈30 × 106 cells) resolution grids were considered. For brevity,
the results presented here only correspond to the medium
resolution grid. Geometries and numerical grids are openly
accessible on JAXA’s Aerodynamics Prediction Challenge(APC) website at https://cfdws.chofu.jaxa.jp/apc/grids/nasacrm/APC1-3/wbh/upacs/fastar/.

The pressure coefficient, CP , contours on the main wing at
AoA = 3.7 and 4.8◦ obtained with FaSTAR (left-bottom) arecompared with the experimental time-averaged PSP by Sugioka et al. (2021) (left-top) in figure 1. The CP distributions
at different spanwise sections of the main wing (available athttps://cfdws.chofu.jaxa.jp/apc/apc3) are also compared withthe experiments performed for JAXA’s APC-3 on the same
flow conditions and geometry of Sugioka et al. (2021) for the
AoA = 4.8◦ . The solutions are in reasonable good agreementwith the experiments and are used as base flows for the GSA.

For the medium resolution grid, a total of 16 AoAs arespanned between 1.5 and 4.8◦ . The RANS solutions (later
used as base flows for the GSA) are first compared with theexperimental time-averaged PSP and pressure taps measurements. The buffet onset is then calculated by using two criteria, namely the lift-coefficient method and the GSA study.

Base flow solutions

3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 3.

GSA spectrum and unstable mode at AoA = 3.8◦ . A comparison with the experimental CP fluctuations is given on the right.

Figure 4.

GSA spectrum and unstable modes at AoA = 4.8◦ . The experimental PSD of CP is given on the right.
Sugioka et al. (2021) found a buffet onset at AoA = 3.7◦ . The
strain gauge buffet onset is selected to be compared with the
GSA based predictions. Similarly to what concluded in theexperiments, GSA predicts buffet onset around AoA= 3.75◦ .
Near the onset at AoA= 3.8◦ , GSA predicts only one unstablemode at f = 615 Hz (St ≈ 0.3), as shown in figure 3. In close
agreement with the experimental CP fluctuations, this unstable mode corresponds to the outboard convection of 3D buffet
cells. The mode shape is also in very good agreement with
the GSA results by Timme (2020). Farther from the onset
at AoA= 4.8◦ , GSA predicts several unstable modes in thesame frequency bump range found in the experiments (figure 4). As well as the buffet cells mode, a second type of highfrequency mode is localized in the wing-junction separation
region. A 2D low-frequency mode at f = 115 Hz (St ≈ 0.06) isLift-coefficient based buffet onset prediction
The lift-coefficient method is a criteria commonly usedin the industry and consists of determining the buffet onset bydetecting a change of the lift-coefficient curve gradient as afunction of the AoA. As shown in figure 2, a linear fit (dashedline) can be used to approximate the lift-coefficient at low angles. A shift of ∆ = 0.1◦ to the right hand side of the linear
fit is then used to determine buffet onset. Both experiments(black empty circles) and RANS solution (red empty circles)
predict buffet onset around AoA = 3.5◦ .

GSA based buffet onset prediction
A second criteria commonly used in the industry is basedon the calculation of a buffet intensity coefficient based on experimental strain gauge measurements. By using this method,

412th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
Table 1. Buffet onset AoA calculated with different criteria
for both experiments (EXP) and numerical simulations (CFD).
Lower and upper limits are reported to account for the sensitivity tests. The experimental results are taken from Sugioka
et al. (2021).

EXP

CFD

Lift

3.50◦

3.50 − 4.00◦

Strain

3.70◦

×

GSA

×

3.75 − 4.15◦

ONSET

3.50 − 3.70◦

3.50 − 4.15◦

CONCLUSIONS
The work here presented is the first fully-3D GSA investigation aimed at characterizing buffet over a large rangeof AoAs and predicting its onset. By calculating the liftcoefficient with JAXA’s RANS solver, the widely used in theindustry lift-coefficient criteria is used to predict buffet onset.
The GSA based buffet prediction is then compared with theexperimental strain gauge measurements carried out at JAXAJTWT1 wind tunnel (Sugioka et al., 2021). An overall good
agreement exists with the experimental findings, although thenumerical agreement with the PSP-based perturbation distributions deteriorates farther away from the onset.
The agreement with the experiments shows the possibility
to use GSA as a computationally cheaper alternative to moreclassical approaches (i. e. URANS) for buffet onset prediction.
To obtain converged information on buffet onset, the estimatedcomputational time saving with GSA with respect to URANS
is in fact of about one order of magnitude. By including sensitivities to grid and turbulence modeling, the GSA-based buffet
onset prediction can be a valuable tool during both design andcertification phases of commercial aircrafts.

Figure 5. Comparison between the experimental CP,rms andthe GSA unstable mode eigen-pressure amplitude at AoA =3.7◦ (top) and 4.8◦ (bottom).

also found. By looking at the experimental root-mean-squarepressure coefficient fluctuations CP,rms in comparison with the
eigen-pressure amplitude distribution associated with the GSA
unstable buffet cells mode (figure 5), it should be however
noted that further away from the buffet onset the agreement
between experiments and GSA deteriorates. Near the onset,
the shock oscillation amplitude is small and the RANS solution used for the GSA study well approximates the timeaveraged flow field. However, for higher angles the shock oscillation amplitude increases and the difference between thesteady RANS solution and the time-averaged flow field result
in larger discrepancies.

ACKNOWLEDGMENTS
The numerical work was partially supported by
JSPS KAKENHI Grant-in-Aid for Early-Career Scientists
20K14953.

Sensitivity tests
REFERENCES

The sensitivity of RANS solutions and GSA results togrid resolution and turbulence model is investigated. For thegrid sensitivity, coarse (≈ 2.8 × 106 cells), medium (≈ 9 × 106
cells) and fine (≈ 30 × 106 cells) grids are investigated. For the
turbulence model sensitivity, the SA-R-QCR200 is comparedagainst the 2-equation Menter Shear Stress Transport (SST)model (Menter, 1994). Additionally for the GSA study, thesensitivity to the perturbation amplitude ε (see section 2.2) isalso analyzed by selecting ε = 1 × 10−1 , ε = 5 × 10−2 andε = 1 × 10−2 . The RANS solutions (and therefore the liftcoefficient based buffet onset prediction criteria) are nearly insensitive to grid resolution and turbulence model tests. For
the GSA study, while turbulence model and perturbation amplitude show no significant sensitivity, the buffet onset is predicted at AoA = 3.75◦ for both coarse and medium resolution
grids and at AoA = 4.15◦ for the fine resolution grid. Table 1
shows a summary of the experimental and numerical onset
predictions using the different methods described above. The
sensitivity tests are accounted for by reporting the onset AoA
lower and upper limits.

Akervik, E., Brandt, L., Henningson, D. S., Hoepffner, J.,
Marxen, O. & Schlatter, P. 2006 Steady solutions of theNavier-Stokes equations by selective frequency damping.
Physics of Fluids 18, 068102.
Arnoldi, W. E. 1951 The principle of minimized iterations inthe solution of the matrix eigenvalue problem. Quarterly of
Applied Mathematics 9.
Bagheri, S., Åkervik, E., Brandt, L. & Henningson, D. S. 2009
Matrix-free methods for the stability and control of boundary layers. AIAA Journal 45.
Barkley, D., Blackburn, H. M. & Sherwin, S. J. 2008 Direct
optimal growth analysis for timesteppers. International fornumerical methods in fluids 57, 1435–1458.
Crouch, J. D., Garbaruk, A., Magidov, D. & Travin, A. 2009
Origin of transonic buffet on aerofoils. Journal of Fluid Mechanics 628, 357–369.
Crouch, J. D., Garbaruk, A. & Strelets, M. 2019 Global instability in the onset of transonic-wing buffet. Journal of Fluid
Mechanics 881, 3–22.

5

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
Paladini, E., Beneddine, S., Dandois, J., Sipp, D. & Robinet, J.-Ch. 2019 Transonic buffet instability: From twodimensional airfoils to three-dimensional swept wings.
Physical Review Fluids .
Richez, F., Lequille, M. & Marquet, O. 2016 Selective frequency damping method for steady rans solutions of turbulent separated flows around an airfoil at stall. Computers
and Fluids 132, 51–61.
Sansica, A., Hashimoto, A., Koike, S. & Kouchi, T. 2022 Sidewall effects on the global stability of swept and unswept
supercritical wings at buffet conditions. , vol. AIAA Paper
2022-1972.
Sansica, A., Robinet, J.-C., Alizard, F. & Goncalvez, E. 2018
Three-dimensional instability of a flow past a sphere: Machevolution of the regular and hopf bifurcations. Journal of
Fluid Mechanics 855, 1088–1115.
Sharov, D. & Nakahashi, K. 1998 Reordering of hybrid unstructured grids for lower-upper symmetric Gauss-Seidel
computations. AIAA Journal 36, 484–486.
Spalart, P. R 2000 Strategies for turbulence modelling and simulation. International Journal of Heat and Fluid Flow 21,
252–263.
Spalart, P. R. & Allmaras, S. R. 1992 A one-equation turbulence model for aerodynamic flows. In 30th Aerospace Sciences Meeting and Exhibit, Aerospace Sciences Meetings, ,vol. AIAA Paper.
Sugioka, Y., Nakakita, K., Koike, S., Nakajima, T., Nonomura,T. & Asai, K. 2021 Characteristic unsteady pressure field ona civil aircraft wing related to the onset of transonic buffet.
Experiments in Fluids (62:20).
Tezuka, A. & Suzuki, K. 2006 Three-dimensional global linear
stability analysis of flow around a spheroid. AIAA Journal
44.
Timme, S 2020 Global instability of wing shock-buffet onset.
Journal of Fluid Mechanics 885.
Visbal, M. R. & Gordnier, R. 2000 A high-order flow solver
for deforming and moving meshes. , vol. AIAA Paper 20002619.

Dacles-Mariani, J, Zilliac, G. G., Chow, J. S. & Bradshaw, P.
1995 Numerical/experimental study of a wingtip vortex in
the near field. AIAA Journal 33, 1561–1568.
Edwards, W. S., Tuckerman, L. S., Friesner, R. A. & Sorensen,
D. 1994 Krylov methods for the incompressible NavierStokes equations. Journal of Computational Physics 110,
82–101.
Giannelis, N. F., Vio, G. A. & Levinski, O. 2017 A review
of recent developments in the understanding of transonicshock buffet. Progress in Aerospace Sciences 92, 39–84.
Guiho, F., Alizard, F. & Robinet, J.-C. 2016 Instabilities in
oblique shock wave/laminar boundary-layer interactions.
Journal of Fluid Mechanics 789, 1–35.
Hashimoto, A., Murakami, K., Aoyama, T., Ishiko, K.,
Hishida, M., Sakashita, M. & Lahur, P. 2012 Toward
the fastest unstructured CFD code ’FaSTAR’. AIAA Paper
2012-1075 .
Ishida, T., Hashimoto, A., Ohmichi, Y., Aoyama, T. &Takekawa, K. 2017 Transonic buffet simulation overNASA-CRM by unsteady-FaSTAR code. AIAA Paper
2017-0494 .
Lee, B. H. K. 2001 Self-sustained shock oscillations on airfoils at transonic speeds. Progress in Aerospace Sciences37, 147–196.
Lehoucq, R. B., Sorensen, D. C. & Yang, C. 1997 ARPACK
user’s guide: solution of large scale eigenvalue problems
with implicitly restarted Arnoldi methods. Technical Note .
Loiseau, J.-C., Robinet, J.-C., Cherubini, S. & Leriche, E.
2014 Investigation of the roughness-induced transition:
global stability analyses and direct numerical simulations.
Journal of Fluid Mechanics 760, 175–211.
Menter, F. R. 1994 Two-equation eddy-viscosity turbulencemodels for engineering applications. AIAA Journal 32,
1598–1605.
Mettot, C., Renac, F. & Sipp, D. 2014 Computation of eigenvalue sensitivity to base flow modifications in a discrete
framework: Application to open-loop control. Journal ofComputational Physics 269, 234–258.

6