Buffet boundary prediction using RANS-based criteria and adjoint methods

Article history: Received 13 January 2022 Received in revised form 20 April 2022 Accepted 24 May 2022 Available online xxxx Communicated by Grigorios Dimitriadis In this work a procedure is proposed for a fast estimate of the transonic buffet boundary on aerofoils and wings. This algorithm aims to be an alternative to unsteady computational fluid dynamics simulations that are prohibitively expensive for complex flow cases. The algorithm combines steady Reynoldsaveraged Navier-Stokes simulations with an adjoint method. Reynolds-averaged Navier-Stokes simulations were used to determine the buffet onset by means of suitable criteria, while the adjoint method was used to compute the sensitivities of the buffet parameters to the flight conditions, and draw the buffet boundary. The procedure was tested for two-dimensional configurations with available experimental data and for two different criteria. Particular attention was paid to the efficiency of the method and its ability to reduce computational cost with respect to purely RANS criteria for buffet boundary estimation. Possible variations to the procedure were also discussed in the results section together with the effect of some key parameters of the proposed algorithm. © 2022 The Author(s). Published by Elsevier Masson SAS. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Transonic buffet is defined as a self-sustained shock oscillation occurring at transonic conditions, and is associated with a dynamic structural response. This oscillation may be detrimental in terms of structural integrity or flight controllability. The exact buffet mechanism is not fully understood. A number of researchers tried to shed light to the mechanism generating and driving the oscillations over aerofoils and wings [1][2][3][4][5]. For an introduction to the topic, the reader is referred to the review paper of Lee [6] and, for recent developments, to the work of Giannelis et al. [7]. The first explanation of buffet onset was given more than fifty years ago by Pearcey [1] who related the oscillation to bubble bursting due to an increase of Mach number or angle of attack. When the shock is strong enough to cause separation extending from the shock foot to the trailing edge (type 3 separation in the characterisation of Mundell and Mabey [8]), low frequency, high amplitude fluctuations are measured in the separated region, together with a periodic and self-sustained shock motion. In recent years, that hypothesis was almost completely surpassed. Nevertheless, the re-lationship between buffet and the aforementioned parameters is not clear, and was studied in many recent works. Tijdeman [2] proposed a characterisation of the shock motion based on flapequipped aerofoils and generalised to plunging or pitching aerofoils. The type A buffet consists of a sinusoidal shock oscillation over the suction side of the aerofoil, with a varying strength and position of the shock throughout the buffet cycle. Type B follows the same pattern of type A, with the shock disappearing during the downstream excursion. Type C sees a strengthening and weakening of the shock as it moves upstream and the propagation into the oncoming flow as a free shock-wave. One of the most plausible explanations behind buffet is the acoustic feedback mechanism proposed by Lee [3] for the shock motion of type A. The shock oscillation on the upper side of the aerofoil generates pressure waves propagating through the separated flow region extending from the shock to the leading edge. Once the leading edge is reached, another disturbance propagates backwards at the local speed of sound. These waves interact with the shock and transfer the energy required to sustain the oscillation. The buffet period is the time necessary for a pressure wave to depart from the shock and reaches the trailing edge plus the time needed for the disturbance to move backwards and hit the shock. In the experimental work of Jacquin et al. [4] on the supercritical aerofoil OAT15A, the authors concluded that there are also upstream travelling perturbations on the pressure side that are diffracted at the leading edge, and play Nomenclature C p pressure coefficient cost function in eq. ( [9] seems to confirm this latest findings. The significant impact of acoustic waves on the shock oscillation over aerofoils was also shown by the more recent experiments of Feldhusen-Hoffmann et al. [10]. Starting from Crouch et al. [5], buffet has been associated with a global instability mechanism, and studied by means of linear stability analyses and the URANS equations. Buffet is seen as a Hopf bifurcation, for which the least stable eigenvalue of the associated linear system crosses the imaginary axis of the complex plane and becomes unstable. The eigenmode associated with the unstable eigenvalue is qualitatively different from the mechanism of Lee [3] and similar to that from Jacquin et al. [4]. A pressure disturbance generates at the shock foot and moves upward along the shock up to the end of the supersonic region. As the perturbation moves upward, the shock approaches the trailing edge and intensifies. A pressure wave generates, goes around the trailing edge, propagates forward along the pressure side and, once at the leading edge, is ingested into the sonic zone. From the same authors, the same approach has been tested on the OAT15A aerofoil [11] and the NACA0012 [12] aerofoil, for which the comparison was done with respect to the early experiments of McDevitt and Okuno [13]. In their study of the flow around the NACA0012 aerofoil, by means of RANS simulations, a strong link between the buffet onset and the appearance of an unstable global mode was found. The same approach was followed by Sartor et al. [14]. These latter also solved the adjoint equation to study the receptivity of the buffeting flow to the location of flow control. The distinction between global modes and adjoint modes emphasizes where the unstable mode is more energetic and where it is more receptive to flow control, respectively. As stated before, buffet is caused by an increase in the angle of attack or the Mach number. These two parameters are used to define the buffet boundaries, upper and lower, where buffet takes place. For a two-dimensional case, Giannelis et al. [15] investigated the effect of both the parameters by means of 2D URANS, while a similar study was carried out by Masini et al. [16] for an aircraft wing. In three-dimensional configurations, the shock dynamics is affected by several factors and presents differences from the twodimensional case. Experimental investigations [17][18][19][20][21] have been carried out on transonic wings, combining steady and unsteady pressure measurements with oil-flow visualisazions and, where possible, particle image velocimetry or laser doppler velocimetry. The experimental findings were corroborated and complemented by subsequent CFD computations. Among them the computational work of Ionovich and Raveh [22] first focused on differences between straight and swept wings and studied the influence of the sweep angle on the shock dynamics. Following the combined effort put into experimental and computational studies, it was found that the 3D buffet is characterised by a broadband spectrum of frequencies rather than a single frequency driving the oscillation. Moreover, the spanwise organisation of the flow in "buffet cells" was shown, characterised by regions of alternated pressure propagating in the spanwise direction, and usually outboard. These features keep constant across different experimental databases and configurations, as shown in the comparative work of Paladini et al. [23]. With the recent advent of pressure sensitive paint (PSP) techniques [21,24], it became possible to acquire snapshots of pressure distribution over the wing and perform modal analysis. In the same way as in some computational works [25,26], Masini et al. [16] employed proper orthogonal decomposition and dynamic mode decomposition to the unsteady PSP measurements of Lawson et al. [21] and showed that the pressure fluctations propagate both inboard and outboard.
Crouch et al. [11] and Timme [27,28] tried to link the buffet onset on 3D wings to a global unstable aerodynamic mode, in the same way as in the 2D case [12]. In the works of Timme [27,28], the onset of buffet on a NASA Common Research Model (CRM) was studied by means of global stability. The base flow was evaluated with a RANS approach, and the Jacobian matrix associated with the spatial discretisation operator was evaluated around the base flow to build the eigenproblem. After the linearisation of the system, the eigenvalues and eigenmodes were analysed to enforce the theory that also 3D buffet can be connected to a global instability mechanism. As the angle of attack increases, the eigenvalues cross the imaginary axis and give rise to the buffet. The study also revealed a buffet cell pattern similar to that found in the numerical simulations of Iovnovich and Raveh [22] and in the experimental work of Dandois et al. [18]. The frequencies on the unstable eigenvalues were distributed around St = 0.4, as shown in the work of Dandois et al. [18]. The stability analysis of Crouch et al. [11] gave light to two classes of modes. The first type of mode, the "oscillating" one, drives the instability as known in the 2D case, while the other one, "stationary", is responsible for adding the 3D component. Their study was performed on an unswept configuration and repeated for a swept wing. The addition of the sweep angle broke the symmetry and the stationary modes became travelling outboard. It is not clear if the 3D buffet can be related to a global instability mechanism.
In parallel, attention was paid to the implications of transonic buffet in flight. Sometimes, as in the case of gusts or emergency manoeuvres, an aircraft can cross the boundary and undergo buffet, and the flight safety can be significantly affected. Although buffet may not compromise the structural integrity of the aircraft, difficulties in piloting may be encountered, and a prompt action of the pilots is required. Therefore, the aircraft must be free from vibration or buffeting under any operating conditions. The buffet boundary can be embodied in the flight envelope so that a threshold velocity is not exceeded, and safety is guaranteed. The (lower) buffet boundary can be defined as a line in the α − M ∞ and C L − M ∞ planes. This coincides with the buffet onset incidence at a given Mach number (or vice versa) and represents the boundary that must not be exceeded in flight. The buffet boundary can also be related to that line separating the regions of attached or partially separated flow, and that of totally separated flow (i.e. from the shock foot to the trailing edge). The constraint on buffet must be accounted for in the aerodynamic optimization process (see. [29]), therefore it arises the need for an accurate method for the buffet boundary determination across a wide range of operating conditions.
In this regard, engineering criteria have been proposed for the buffet boundary estimation, based on the results of experimental investigations. Because of the extension of the separated region from the shock foot to the trailing edge, a common indicator based on trailing edge measurements, is the trailing edge pressure divergence [1]. An alternative quantity, also based on pressure measurement over the aerofoil, is the unsteady normal force, whose divergence corresponds to the buffet onset. Alternative criteria mentioned in the paper of Mabey [30] consist of methods using a wing-root strain gauge or methods evaluating the breaks in the C L − α and C D − C 2 L curves. These, however, compared to an actual buffet boundary, did not prove to be good indicators of buffet. Lawson et al. [21] assessed the performance of the aforementioned criteria, together with two other, that were based on tip acceleration and pitch moment break, using experiments conduced by ARA for the RBC12 transonic wing. They underlined that different criteria may provide quite different predictions. Using computations, Chung et al. [31] employed some of the most common criteria in conjunction with steady and unsteady RANS and compared the results with the experiments on a NACA0012 aerofoil from McDevitt and Okuno [13]. Although not far from measurements, all criteria overestimated the onset angle of attack of the buffet, and the accuracy of some reduced as the Mach number increased. Unsteady simulations gave better results, but with a significant increase in computational costs. Since the proposed criteria require the computation of the steady flow at several angles of attack for each Mach number, the overall procedure may not be affordable when coming to three-dimensional cases. A number of authors employed URANS to predict the buffet for both two dimensional [32][33][34][35][36][37] and three-dimensional cases [22,38]. Being unsteady RANS simulations influenced by many factors, such as turbulence modelling, numerical schemes, spatio-temporal discretisation and influence of the wind tunnel geometry, no consensus was found on their ability to adequately describe the buffet phenomenon and predict the onset. Therefore, an increasing number of authors have started adopting hybrid RANS/LES approaches to analyse the flow around two-dimensional [25,33,39,40] and three-dimensional con-figurations [41][42][43][44], or even LES simulations [9,45]. Although they showed significant improvements in the characterisation of buffet, scale-resolving simulations are too expensive to be applied to complex three-dimensional configurations or to be repeated at different conditions. This raises the need for a reliable but affordable method to predict the buffet onset that can be repeated at several flight conditions.
In this paper, the concept of RANS-based criteria is exploited to formulate an algorithm for a fast buffet boundary evaluation. For this purpose, RANS simulations are performed in conjunction with an adjoint method implemented in the Helicopter Multi-Block (HMB3) solver. Section 2 is devoted to the description of the computational apparatus. Section. 3 shows the results obtained with the procedure and the comparison with purely RANS-based criteria, with particular attention to the some parameters involved, variations of the algorithm and the relative computational times, before drawing some conclusions in section 4.

Computational model for fluid flow
Numerical simulations of the steady flow fields have been performed using the Helicopter Multi-Block (HMB3) [46,47] flow solver, a three-dimensional, fully implicit, structured, multi-block code for the solution of the Navier-Stokes equations. The Navier-Stokes equations are discretised using a cell-centered finite volume approach. The computational domain is divided into a finite number of non-overlapping control volumes, and the governing equations are applied to each cell in turn. Also, the Navier-Stokes equations are re-written in a curvilinear co-ordinate system which simplifies the formulation of the discretised terms since bodyconforming grids are adopted here. The spatial discretisation of equation leads to a set of ordinary differential equations in time, where W and R are the vectors of cell conserved variables and residuals respectively and V is the cell volume. The convective terms are discretised using Osher's upwind scheme [48]. A monotone upstream-centered scheme for conservation laws (MUSCL) variable extrapolation [49] is used to provide second-order accuracy with the Van Albada limiter [50] to prevent spurious oscillations around shock waves. The solver offers several one-, two-, three-, and four-equation turbulence models. In addition, LES, DES, delayed DES (DDES), improved DDES (IDDES), SAS and PANS methods are also available.
The integration in time of eq. (1) to a steady-state solution is performed using an implicit time-marching scheme by where n + 1 denotes the time (n + 1) * t. Eq. (2) represents a system of non-linear algebraic equations and to simplify the solution procedure, the flux residual R ijk W n+1 ijk is linearised in time as follows, where (2) now becomes the following linear system The linearised system of equations is solved using the generalised conjugate gradient (GCG) method with a block incomplete lowerupper (BILU) factorisation as a pre-conditioner [51]. The Jacobian is approximated by evaluating the derivatives of the residuals with a first-order scheme for the inviscid fluxes. The first-order Jacobian requires less storage and ensures a better convergence rate to the GCG iterations. The steady-state solver for turbulent flows is formulated and solved in an identical manner to that of the mean flow. The eddy-viscosity is calculated from the latest values of the turbulent variables, e.g. k and ω, and is used to advance the mean and the turbulent flow solutions. An approximate Jacobian is used for the source term of the models by only taking into account the contribution of their dissipation terms, i.e. no account of the production terms is taken on the left-hand side of the system.

Adjoint formulation
To compute the derivatives of global, aerodynamic coefficients with respect to flight parameters, the sensitivity equation cast in adjoint form is solved. The underlying idea is to write explicitly the cost functional I as a function of the flow variables (W ) and of the design variables α, that is, I = I(W(α), α). The flow variables are subject to satisfy the fluid dynamics governing equations written in compact form as Formally, taking the derivative of I with respect to α we obtain: By introducing the adjoint variable λ as the solution of the following linear system: equation (6) can be rewritten as: To obtain eq. (8) the governing equations (5) have been differentiated: The computational cost of the dual sensitivity problem (7)- (8) scales with the number of outputs, since the right-hand side of (7) depends on I , but it is independent of the input parameters. The adjoint form of the sensitivity equation is particularly efficient for cases where the number of (output) cost functionals is small, while the number of (input) design variables is large. Options for the solution of the linear system in eq. (7) are a fixed-point iteration scheme or a nested Krylov-subspace method based on the Generalized Minimum Residual (GMRES) method. For further details the reader is referred to the work of Biava and Barakos [52] and Biava et al. [53].

RANS-based criteria
The RANS-based criteria considered in this study are now described. Similar correlations can be found in the works of Chung et al. [31] and Lawson et al. [21]. The first indicator can be found looking at the change in the lift coefficient slope, as suggested by Pearcey and Holder [54]. This relates to the effect of the flow separation at the shock foot when the shock moves downstream as the angle of attack increases. To determine the buffet onset from the lift curve, the straight line approximating the linear regime of the lift curve is shifted by α = 0.1 deg and the buffet onset is identified as the intersection between the line and the lift curve itself (see Fig. 1, left).
The second criterion is based on the slope of the moment coefficient. The pitching moment variation is sensitive to the length between the reference point and the location of the separated region. This length acts like an amplifier of the changes in the mean aerodynamic loads. Indeed, as the separated region moves aft with the movement of the shock, this effect is accentuated. The deviation point from the linear range is seen as buffet onset [31]. An example is shown in Fig. 1 (centre).
The same procedure is used for the trailing edge pressure divergence, as shown in Fig. 1 (right). For transonic aerofoils, the pressure coefficient changes as a result of the flow separation due to the shock on the suction side. This coincides with a lift drop and the onset of buffet. As specified by Chung et al. [31], there is no value for the exact deviation from the linear regime, but an offset of C p = 0.05 was recommended in Ref. [21,31].
The accuracy of the presented criteria is limited by the performance of RANS simulations that in the past have shown different results according to the turbulence model, numerical scheme and spatio-temporal discretisation adopted.

Algorithm formulation
In this section, the algorithm used in this work is presented. The present algorithm makes use of the adjoint method implemented in the solver HMB3 to reduce the CPU cost associated with the aforementioned RANS-based criteria and was not intended to improve the accuracy in the prediction of the buffet boundary. Indeed, the accuracy of the algorithm is strongly related to that of the criteria described in sec. 2.3. A flow-chart of the algorithm, implemented in MATLAB, is shown in Fig. 2 In the current study of the buffet boundary, we will account only for Mach number and angle of attack changes (an alternative might be the lift coefficient). A third parameter like the Reynolds number can also be accounted for but it would increase the required CPU time.
The Mach number is an important parameter in defining the flight conditions and we opted to keep it fixed and evaluate the buffet onset in terms of the angle of attack. It must be pointed out that the opposite is also possible. The range of flight conditions we are interested in is represented by an interval of Mach numbers The selected criterion must be as accurate as possible over the entire interval of angles of attack. The efficient procedure employed here uses RANS-based criteria. A generic buffet parameter C B must be introduced. In the case of the aforementioned criteria, some possible candidates are the aerodynamic coefficients (C L or C M ), the trailing edge pressure coefficient, or one of their derivatives with respect to the flight parameters. Then, a target value C * B for the buffet onset parameter C B must be specified as a function of the Mach number. C * B is defined as the value assumed by C B at the onset. The choice of a suitable buffet parameter is critical and will determine the outcome of the procedure. Moreover, the buffet onset parameter must be inferred from steady computations and, therefore, dynamic parameters must be avoided. The adjoint method used to move along the buffet boundary can also be used in conjunction with unsteady simulations. This is possible using parameters like the shock oscillation frequency or amplitude, and the higher order moments of the flow variables. Nevertheless, the present method is meant to be the least expensive possible, and the use of unsteady simulations increases the required CPU time. Once an increment in the Mach number M ∞ is fixed, the point is propagated along the buffet boundary where C B follows: The next point is evaluated increasing (or decreasing) the Mach number by the fixed increment M ∞ and by updating the angle of attack according to eq. (10): where C B,α and C B,M ∞ were now used to express the partial derivatives in eq. (10). iThe angle of attack is then corrected to reach the target value of C B through: To ensure that, in the correction, the sensitivities do not change significantly, a tolerance α is introduced. If the correction exceeds the range determined by α , the sensitivities are re-computed and a new correction is made; otherwise, the Mach number is updated and the next point is chased. The process is initiated at the two extremes of the Mach interval and the buffet boundary is traced going towards the centre of the interval. If more than one corrections are needed, a higher order expansion can be used to speed up the process and increase the accuracy of the procedure. In that case, when two corrections are applied, i.e. C B , C B,α are known at two different angles of attack, eq. (12) can be replaced by a quadratic expansion, where the second-order derivative is assumed constant between the two known points. The same can be done for higher order expansions. Alternatively, when the number of sample points is high, the C B − α curve can be interpolated through a high-order polynomial without relying on the result of the adjoint method (since the derivatives are no longer needed). In this case, a purely RANS-based criterion is used to estimate the buffet onset. All these alternatives will be explored in section 3 where results of the proposed method are presented.

Test cases description OAT15A Wing Section
The first configuration is the supercritical aerofoil OAT15A, studied in an experimental campaign by Jacquin et al. [4,55]  For the case at M ∞ = 0.73 and α = 3.5 deg, a laser Doppler velocimetry (LDV) system was used to acquire velocity-field data and compute statistics. At all other flow conditions the pressure measurements were complemented with oil-flow and schlieren visualisations.

NACA0012 Wing Section
The experiments used for the flow around the NACA0012 section were performed in the Ames High Reynolds Number Facility [13]. Although these experiments are more than 30 years old, they still represent the broadest database about buffet onset, covering a wide range of flight conditions. In particular, Mach numbers are in the range of 0.71-0.8 and Reynolds number span between 1 and 10 million. The measurements of both steady and unsteady pressure measurements, allowed to detect the buffet onset at several conditions. The tunnel walls were adapted to follow the free air streamlines, while the sidewall interference was reduced by thinning the sidewall boundary layer by means of suction applied on porous panels. The experiments were supported by a shadowgraph system. A more detailed description is given in reference [13].

RANS computations and comparison with experiments
The computational grids considered have a typical C-H topology. The 2D grid consists of 320 or 360 cells around the aerofoil, Update M i+1   a more upstream position of the shock due to the presence of the sidewalls. Since the exact prediction of the shock position is beyond the scope of this work, the current agreement is seen as satisfactory.

Adjoint method validation
The NACA0012 at Re = 6.0 × 10 6 was taken as the reference case for comparing the performance of the adjoint methods in evaluating sensitivities with respect to the flight conditions. Fig. 6 shows the comparison between the derivative ∂C L ∂α evaluated by means of the adjoint method (denoted as ADJ) and finite differences (denoted as FD). The overall agreement is good at every Mach number. Some small differences arise when the angle of attack is far beyond the experimental buffet onset (black vertical lines in Fig. 6) and some oscillation are present, e.g. at M = 0.77, 0.79. This is due to limitations with finite differences and the employed steady flow approximation. Overall, we do not expect to cross this threshold, and the agreement is, overall, good. The same can be stated for the derivative Fig. 7. Again the region where values differ is well beyond the buffet onset, the flow in unsteady, and will not be used in the buffet onset estimation.

Choice of the onset criterion
Among the criteria mentioned in the previous sections, the ones based on C L and C M were chosen to test the algorithm. The main reason is that integral aerodynamic coefficients can be efficiently used in conjunction with the adjoint method. Indeed, the criterion based on the trailing edge pressure coefficient was not used here. The adjoint variable in eq. (7) is introduced to avoid the evaluation of the derivatives of the flow variable with respect to the flight parameters. Building the algorithm with ∂C p,T E /∂α as the buffet parameter, as an example, would reintroduce the issue of evaluating such terms and the only possible way would be to rely on finite differences. Conversely, assuming C B = ∂C L /∂α or C B = ∂C M /∂α allows us to use eq. (6) to compute the derivatives.  Moreover, the target buffet indicator must be easy to estimate all over the interval of interest. It must also show some trend that is independent on the geometry under analysis. Fig. 8 shows the behaviour of the two derivatives at the buffet onset predicted by the two criteria for two different geometries and Reynolds numbers. The target buffet indicator can be, in both cases, approx-imated with a linear function once the values at the minimum and maximum Mach number are known. We see that this approximation is more correct for C B = ∂C L /∂α (solid lines in Fig. 8  buffet boundary in comparison with the experiments, as shown in Fig. 9.

Results and discussion
Here the results obtained for the NACA0012 at Re c = 6 ×10 6 are analysed in detail. Fig. 10 shows the comparison between different runs corresponding to those listed in Table 1 When the tolerance α is decreased, the algorithm takes more time to converge like in the switch from case (a) to case (b).
Moreover, using too fine tolerance on α, possibly emphasizes the differences between the approximated target buffet coefficient and the actual value, hence there is no guarantee that decreasing the tolerance under a certain threshold brings significant benefits in the estimate. Increasing M ∞ , the role of the adjoint in the evaluation of the derivatives becomes more important. If M ∞ is small, an error made in the approximation is compensated by the linear correction of eq. (12), as long as the distance between the buffet onset and the current point α i is relatively small. Increasing the Mach number increment, the error propagates and the first investigated point may be too far from the actual buffet onset. If this happens, the linear extrapolation may take many iterations or reach regions where the adjoint method no longer provides good estimates of the derivatives.  To overcome this drawback, some modifications were introduced. To avoid to fall in regions where the derivative computed with the adjoint method has wrong value and sign, an angle restriction (conf. in Table 1) was applied. For conventional aerofoils, and in this range of Mach numbers, the buffet onset angle of attack is lower with higher Mach numbers, and it is, therefore, reasonable to estimate the buffet boundary as the straight line passing through the extreme points. From that estimate, the angle of attack is confined with a certain interval, and 1-2 degrees are a good compromise to avoid falling away from the onset without preventing non linearities of the criterion. The effect of the angle confinement can be seen comparing cases (d) and (e). In case (e), the number of additional runs required is decreased. Moreover, if the number of attempts increases, a quadratic or higher order interpolation can be used to reduce the number of points required. From two points, given the value of C B and its derivatives, a second order interpolation can be done and the order increases with the number of attempts. Nevertheless, using a too high order of the interpolation may generate many different roots in the interval of interest, so in this work a maximum order of 4 was used.  Using a higher order interpolation, the number of runs required decreased. The procedure in cases (d-f) required a smaller number of runs with respect to the linear case (c). Also, when the number of attempts increases, an alternative to interpolating C B is to reconstruct directly C L (α) or C M (α). While the points used to interpolate C B are affected by the errors introduced by the adjoint evaluation, the values of the aerodynamic coefficients are in principle more accurate, coming from the same RANS computations used to evaluate the onset in the basic version of the criteria. Using this approach (C L -based or C M -based in Table 1), the number of simulations required is further reduced and the overall accuracy is increased. From now on, the results presented were obtained employing this last approach.
In this case (f), the number of simulations required to draw a buffet boundary is lowered by around 60%, at the cost of a mean error of 0.074 deg and maximum error of 0.256 deg, with respect to the purely RANS-based criterion. Increasing the tolerance α (case (g)), the accuracy of the algorithm is slightly reduced, as well as the CPU cost. It must be pointed out that such an estimate in the reduction of CPU costs is meaningful only if the Mach number interval and the number of intermediate Mach numbers are the same. This also puts into evidence the advantage of the proposed algorithm over the purely RANS-based criteria. While the number of computations required by the purely RANS-based criteria is proportional to the number of intermediate Mach numbers, the computational cost of the algorithm is not so heavily affected by the number of intermediate points. Therefore, the smaller M ∞ , the higher is the saved CPU cost.
The same study is repeated for the OAT15A test case (Fig. 11). The Mach number span from 0.69 to 0.77. Once again, decreasing the Mach number increment does not necessarily result in higher costs, as long as the tolerance on the angle of attack is also increased. In terms of performance, the algorithm introduces a mean error of around 0.1 deg (that is comparable with the toler-ance used) and a maximum error of 0.24-0.31 deg. The number of simulations required is around half of those needed by the purely RANS-based criterion. Here, a slightly higher discrepancy is found with respect to the previous test case. The main reason is that the function C L,α (M) (Fig. 8) slightly departs from the linear behaviour assumed, making the computation less accurate.
The procedure was then applied in conjunction with the C M,αbased criterion for both configurations (Figs. 12 and 13). For case (b), both figures underline the drastic increase in the number of computations required when decreasing the Mach number increment M ∞ , keeping a small tolerance on the angle of attack. Although the accuracy of the computation does not vary between cases (a) and (b), for cases (b) the CPU cost was higher by 70% for the NACA0012 and 27% for the OAT15A.
Overall, this second criterion is slightly more inaccurate, apparently due to the magnitude |C M,α | being smaller than |C L,α |.
The difference in the linear approximation of these quantities with the Mach number is hence more significant in the first case. This results in a slightly higher inaccuracy. On the other hand, this criterion was shown to be closer to the experimental results. Also in this case, the developed algorithm guarantees a reduction of 50-55% of the number of computation required by the purely RANSbased criterion.
It must be pointed out that the number of computations required included the computations needed to evaluate the buffet boundary at the minimum and maximum Mach numbers. In particular, this amounts to 47 computations for the NACA0012 case and 46 for the OAT15A case. This means that the actual number of computations to complete the buffet boundary with the proposed method is actually even smaller.

Conclusions and future work
In this work, an algorithm for a fast prediction of the buffet boundary on aerofoils and wings was developed using steady RANS-based criteria and adjoint methods to compute flow sensitivities. The derivatives C L,α and C M,α were chosen as buffet indicators. Once determined, they can be approximated as linear functions of the Mach number. The algorithm fixes the Mach number and evaluates the onset angle of attack corresponding to the target value for each buffet indicator. The adjoint method was used to evaluate the derivatives with respect to the flight conditions. Other buffet onset criteria, mentioned in sections 1 and 2, were discarded because they involve local quantities, like the trailing edge pressure, that are not as easily available as the integrated loads.
The RANS, 2D computations and the adjoint method results were assessed separately. The RANS results were sufficiently close to experiments. Flow sensitivities calculated with the adjoint method were close to those obtained with finite differences. Some discrepancies were present when the angle of attack was beyond the buffet onset, due to the steady flow approximation used. For that reason, an angle confinement was introduced in the formulation.
The algorithm was tested using two configurations with available experimental results for the buffet onset. The procedure was repeated for both criteria and configurations. The use of linear extrapolation to find the angle of attack was efficient only when small Mach number increments were used. Otherwise, results could depart from the buffet boundary and the costs of the algorithm was increased. High order interpolations gave benefits both in terms of accuracy and cost reduction. Overall, the algorithm allowed to recover the results obtained with purely RANS-based criteria. The algorithm was not formulated to improve the accuracy in the prediction of the buffet boundary with respect to the experiments, but only to improve the efficiency. Indeed, the procedure reduced the computational costs by 50-60% with respect to the aforementioned criteria.
Future work will be devoted to explore the ability of the algorithm to work on more complex geometries e.g., transonic wings. Moreover, the focus will be on modifying the method to use unsteady simulations. Unsteady parameters may give improved predictions of the buffet boundary. Finally, the present study may be enriched with the introduction of a third parameter like the Reynolds number, to widen the range of flight conditions considered.

Declaration of competing interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Andrea Petrocchi was supported by the European Union.