## Abstract

Most forecasts predict an annual airline traffic growth rate between 4.5 and 5% in the foreseeable future. To sustain that growth, the environmental impact of aircraft cannot be ignored. Future aircraft must have much better fuel economy, dramatically less greenhouse gas emissions and noise, in addition to better performance. Many technical breakthroughs must take place to achieve the aggressive environmental goals set up by governments in North America and Europe. One of these breakthroughs will be physics-based, highly accurate and efficient computational fluid dynamics and aeroacoustics tools capable of predicting complex flows over the entire flight envelope and through an aircraft engine, and computing aircraft noise. Some of these flows are dominated by unsteady vortices of disparate scales, often highly turbulent, and they call for higher-order methods. As these tools will be integral components of a multi-disciplinary optimization environment, they must be efficient to impact design. Ultimately, the accuracy, efficiency, robustness, scalability and geometric flexibility will determine which methods will be adopted in the design process. This article explores these aspects and identifies pacing items.

## 1. Introduction

According to Boeing (http://www.boeing.com/boeing/commercial/cmo/) and Airbus (http://www.airbus.com/company/market/forecast/) forecasts, worldwide airline traffic in terms of revenue passenger kilometres is expected to grow 4.7–5% per year in the next 20 years, essentially doubling every 15 years. This rate of growth is much larger than the world gross domestic product growth rate over the same period, forecasted to be around 3.2%. Accommodating such growth will require many new aircraft and more flights, plus new runways, airports and traffic control systems. If the commercial aviation industry remains at the present technology level, the likely consequences are a significant escalation in harmful greenhouse gas (GHG) emissions, unprecedented traffic jams in the sky, non-stop noise near airports and more citizens unhappy about the effects on their health and their quality of life. The consequence of burning fossil fuels is well established in their long-term impact on climate and global warming due to GHG emissions, primary being CO_{2} and NO_{x}. Two recent studies in both the USA and the UK appear to directly link aircraft noise with cardiovascular and other diseases [1,2].

The importance of reducing the environmental impact of aviation has been fully realized by many governments and private industry; and many researchers have started efforts to address the environmental concerns [3–5]. Words such as *green* or *sustainable* are coined to describe future aircraft or aviation with reduced environmental impact. In 2010, the US National Science and Technology Council released the National Aeronautics Research and Development Plan (http://www.whitehouse.gov/sites/default/files/microsites/ostp/aero-rdplan-2010.pdf), which sets the following goals for aircraft in 10 years:

— increase lift/drag ratio by 25%;

— reduce fuel burn by 70% compared with Boeing 737/CFM56;

— reduce noise by 62 dB cumulative below current FAA standard for large subsonic jet aircraft; and

— reduce NO

_{x}emissions by 80% below current international standard.

Many breakthroughs must take place in the coming decades to satisfy the ever more stringent environmental regulations. One of the breakthroughs will be physics-based highly accurate/efficient and robust aircraft and engine design tools, and noise prediction tools. The most critical among them are computational fluid dynamics (CFD) tools capable of handling the entire flight envelope from take-off to landing, and predicting the highly unsteady and turbulent flow inside an engine. At present, most CFD design tools are based on the second-order finite volume method on hybrid unstructured meshes capable of handling complex geometries [6]. The governing equations are the Reynolds-averaged Navier–Stokes equations using a turbulence model such as the Spalart–Allmaras model [7] or detached eddy simulation [8] to handle turbulent flows at high Reynolds numbers. These tools have proved to be very useful in predicting flow at the cruise condition and were used heavily in the design of the latest Boeing and Airbus commercial aircraft. However, they have generally failed to predict highly separated flow for high-lift configurations during take-off and landing, because a statistically steady mean flow may not exist at such flow regimes. In addition, the highly separated turbulent flow is dominated by unsteady vortices of disparate scales, whose accurate resolution calls for high-order CFD methods, at least third-order accurate [9]. Shown in figure 1 is the vorticity distribution computed with a second-order and a fourth-order discontinuous Galerkin (DG) type scheme for an isentropic vortex propagating through the domain eight times. The second-order simulation was performed on a fine mesh, while the fourth-order scheme was conducted on a coarse mesh so that both have the same number of degrees of freedom (*n*_{df}) and roughly the same computational cost. Note that the fourth-order scheme preserves the vortex strength much more accurately than the second-order scheme.

Another major design parameter is aircraft noise, which is composed of three major sources: airframe noise, engine noise and landing gear noise (which is singled out because it is the dominant noise during take-off and landing). It is mostly produced by unsteady, turbulent flow through the engine and around major airframe components. The direct computation of aircraft noise is extremely challenging, as the noise (or pressure fluctuation) is governed by the same equations governing the flow, the Navier–Stokes equations. However, the pressure fluctuations are normally orders of magnitude smaller than the mean flow quantities. To overcome the difficulty, various acoustic analogies [10,11] are used in which the noise computation is decoupled from the flow simulation. From the unsteady aerodynamic flow, acoustic sources are identified, and then propagated to the far field by solving a simpler set of governing equations such as the Euler equations or a linear wave equation. This approach has achieved remarkable success in many applications. Refer to the review article [12] in this Theme Issue for details. The fidelity of the approach hinges upon the following questions:

(1) How accurate is the near-field flow simulation?

(2) Can the noise problem be partitioned into a source and propagation problem?

(3) How accurate can the noise sources be propagated to the far field?

It has been shown in many computational aeroacoustics (CAA) studies that high-order methods with very low dissipation and dispersion errors are critical in noise propagation problems. In addition, many prior noise simulations demonstrated that an accurate near-field computation is often the determining factor of how reliable the noise prediction is. In order to accurately compute noise for the entire flight envelope, accurate near-field aerodynamic simulations call for high-order methods.

## 2. Review of high-order methods

High-order methods have received considerable attention from the CFD community in the past two decades owing to their potential of delivering higher accuracy with lower cost than low-order methods. Before proceeding any further, let us first clarify what we mean by ‘high order’. A numerical method is said to be *k*th-order accurate (or of order *k*) if the solution error *e* is proportional to the mesh size *h* to the power *k*, i.e. *e*∝*h*^{k}. In the aerospace community, high-order means third- or higher-order accuracy.

Many types of high-order methods have been developed to deal with a diverse range of problems. Interested readers can refer to several review articles on high-order methods [13–15]. At the extremes of the accuracy spectrum, one finds the spectral method [16] as the most accurate, and first-order methods (the Godunov method [17], for example) as the least accurate. Many high-order methods have been used successfully in CAA and CFD applications. An incomplete list includes the compact methods [18,19], essentially non-oscillatory (ENO)/weighted ENO methods [20–22], dispersion-relation-preserving (DRP) [23] and upwind DRP methods [24], weighted compact nonlinear schemes [25], spectral element [26], SUPG (streamline upwind Petrov–Galerkin) [27,28], staggered grid multi-domain [29], high-order *k*-exact finite volume [30,31], residual distribution [32], DG [33–36], spectral volume/difference [37–39], PnPm [40] and correction procedure via reconstruction (CPR) [41–43].

Which high-order methods will be embraced by future design tools? Distinguishing features may include geometric flexibility, accuracy/efficiency, robustness and scalability. Note that we put accuracy/efficiency as a single criterion to emphasize that they must be considered together. When we determine which method is the most accurate, we need to look at the solution error produced with roughly the same amount of central processing unit (CPU) time. In order to limit the scope of this paper, we will focus on discontinuous high-order methods capable of handling unstructured meshes.

The ancestor of discontinuous methods was the Godunov finite volume method, in which the numerical solution is discontinuous across element interfaces. In order to ensure conservation, a unique flux must be used at an interface. This flux was computed analytically in the original Godunov method by solving a Riemann problem, but was later replaced with various approximate Riemann solvers [44,45] to reduce the computational cost. The Godunov method is compact in that the scheme residual of an element depends on itself and its immediate neighbours. Compact methods incur minimum amount of data communication on modern parallel computers including CPU and graphics processing unit (GPU) clusters. On a GPU card with thousands of compute cores, it was demonstrated that minimizing communication was paramount in achieving the best performance [46].

Unfortunately, the Godunov method is not accurate at all, being only first-order. The only way to rescue it was to extend it to higher-order accuracy. The most obvious idea was borrowed from a finite-difference method by using neighbouring elements to build a higher-order solution polynomial, as shown in figure 2*a*, resulting in higher-order methods such as the MUSCL [47] and ENO/WENO methods [20]. Later this idea was successfully extended to unstructured meshes and often called *k*-exact finite volume method [30]. The biggest drawback of these methods is that they are not compact. For example, the scheme stencil for a second-order finite volume method includes neighbours' neighbours. The other idea to extend the Godunov method to higher order was borrowed from a finite-element method by defining multiple degrees of freedom (d.f.) on an element, as shown in figure 2*b*. The reconstructed solution polynomial uses only local data from the element itself. The scheme stencil includes itself and only its immediate neighbours. This was exactly how the DG method [35,36] was conceived.

## 3. Recent progress in compact discontinuous methods

To present the basic idea, we consider the following one-dimensional conservation law
3.1
where *Q* is the state variable and *F* is the flux. The computational domain [*a*,*b*] is discretized into *N* elements, with the *i*th element defined by *V*_{i}≡[*x*_{i−1/2},*x*_{i+1/2}]. Each element can be transformed into the standard element [−1, 1] using a linear transformation if necessary. A degree *p* polynomial basis is denoted by .

### (a) Discontinuous Galerkin method

The approximate numerical solution is piecewise continuous with discontinuities at element interfaces. On element *V*_{i}, the solution is approximated with a polynomial of degree *p*, i.e.
3.2
where *u*_{i,j} is the expansion coefficient. Let *W* be a smooth weighting function. Substituting (3.2) into a weighted residual form of (3.1) on *V*_{i}, we obtain
3.3
In order to achieve conservation, the flux across an element interface must be unique. Therefore, the flux in (3.3) is replaced with a unique Riemann flux
3.4
where is the solution on the left side of interface *i*+1/2, and is the solution on the right side of the interface. Equation (3.3) then becomes
3.5
Let W be each of the basis functions . We then obtain just enough equations to update the d.f. As the flux is usually a nonlinear function of the state variable, a Gauss quadrature formula is used to compute the volume integral. Finally denoting *Q*_{h} the global d.f., the equations can be cast in the following form:
3.6
where *M* is the global mass matrix, and *R*_{h}(*Q*_{h}) the global residual.

### (b) Related formulations

The Runge–Kutta DG method was first published in the late 1980s [36]. During the 1990s, the CFD community started to pay attention to the idea. However, the *k*-exact finite volume (FV) method received much more research effort, as unstructured grid methods became the state of the art. Since the 1990s, many researchers have tried to improve the DG method in whatever way possible. For example, a quadrature-free DG approach was developed in [48], and nodal DG (NDG) method developed in [49]. The present author developed an FV version of the DG method and called it the spectral volume (SV) method [15]. The SV method, however, ran into stability and efficiency issues in three dimensions. The search for a more efficient finite-difference-like DG method led Liu *et al*. [37] to develop the so-called spectral difference method. In one dimension, the method is equivalent to the staggered grid multi-domain method [29]. Here is the basic idea.

Given a degree *p* solution polynomial *Q*_{i}(*x*,*t*) on element *V*_{i}, construct a flux polynomial that is one degree higher at *p*+1. Then the solution is updated using an FD-like formula, i.e.
3.7
Two sets of grid points, i.e. the *solution points* and *flux points*, are defined in each element. The flux points must always include the two endpoints so that a unique flux can be imposed at the cell interface to achieve conservation. Let the position vector of the *j*th solution point at cell *i* be denoted by *x*_{i,j}, and the *k*th flux point at cell *i* be denoted by *x*_{i,k}. Denote *Q*_{i,j} the solution at *x*_{i,j}. Given the solutions at the solution points, an element-wise degree *p* polynomial can be constructed using a Lagrange-type polynomial basis, i.e.
3.8
where *L*_{i,j}(*x*) is the cardinal basis function. Next, compute the flux values at the flux points. For the interior flux points, the value is computed based on the state variable at the flux point, i.e. . At the two endpoints, Riemann fluxes are again used. Then the degree *p*+1 flux polynomial is built from the flux values at the flux points
3.9
The SD method ran into stability issues on a triangular mesh [39]. Later, the introduction of a new basis function for the flux appears to have fixed the problem [50].

The search for an FD-like formulation continued and led to the development of the flux reconstruction or CPR method [41]. The difference between the SD and CPR methods lies in how the flux polynomial is constructed. In CPR, the flux polynomial is written as
3.10
where is a degree *p* internal flux polynomial approximating *F*(*Q*_{i}(*x*)) in some sense, and *σ*_{i}(*x*) is the correction flux polynomial of degree *p*+1. One way to compute is
3.11
At both endpoints of element *i*, the flux polynomial should take the values of the Riemann flux, i.e.
3.12
The correction flux polynomial is then further expressed in the following form to satisfy the two boundary conditions:
3.13
where both *g*_{L}(*x*) and *g*_{R}(*x*) are degree *p*+1 polynomials called correction functions, which satisfy
3.14
Equation (3.7) then becomes
3.15
Many correction functions were presented in [41], corresponding to different numerical schemes, including the DG and SD/SV, and other new methods were discovered for the first time. Note that equation (3.15) can be written as
3.16
where is the correction polynomial of degree *p*. Since the pioneering work by Huynh [41], the CPR formulation in the form of (3.16) was first extended to triangular and mixed elements [43,51,52]. Energy stability was proved in [53], and other members of the family were discovered in [54,55]. See [42] for a comprehensive review of recent developments.

The relative performance of quadrature-based DG (QDG), NDG, SD and CPR formulations were compared in [56] by solving a vortex propagation problem governed by the Euler equations on a quadrilateral mesh. Figure 3*a* shows the density error versus work unit, which is a scaled CPU time used in the International Workshop on High-Order CFD Methods [9]. In this idealized comparison, CPR appears to perform the best. Figure 3*b* displays the cost per d.f. for these schemes at different orders of accuracy. It is interesting to note that the cost per d.f. decreases for the CPR scheme with increasing *p*. This is due to the fact that these schemes are one-dimensional in each coordinate direction on a quadrilateral element.

In figure 4, the performance of the second-order finite volume method (*p*=1) in the TAU code was compared with CPR schemes of various orders with a benchmark problem of flow over a bump in a channel, also from the International Workshop on High-Order CFD Methods. The entropy error was plotted against the work units. Note that all high-order CPR schemes (*p*>1) outperformed the second-order FV scheme for this problem.

## 4. Towards a high-order CFD design tool

After decades of research and development mostly in academia and government laboratories, adaptive high-order methods started to attract more attention from industry in the past decade. In the USA, Boeing engineers studied both the SUPG and DG methods [28] for aerodynamic problems. It was concluded that, though these methods demonstrated a lot of potential, much remains to be done for the high-order methods to be used routinely in a design tool. In Europe, the Adaptive Higher-Order Variational Methods for Aerodynamic Applications in Industry (ADIGMA) project [57] supported a consortium consisting of 22 organizations, which included the main European aircraft manufacturers, the major European research establishments and several universities, all with well proved expertise in CFD. The goal of ADIGMA was the development and utilization of innovative adaptive higher-order methods for the compressible flow equations enabling reliable, mesh-independent numerical solutions for large-scale aerodynamic applications in the aircraft industry. ADIGMA's follow-on project Industrialisation of High-Order Methods (IDIHOM) focused on industry problems. In addition, two international workshops on high-order CFD methods were held in 2012 and 2013, respectively. One of the objectives is to provide an open and impartial forum for evaluating the status of high-order methods in solving a wide range of flow problems. Some of the findings are documented in a review paper [9]. In the following subsections, we describe recent progress in several major pacing items, which will determine how soon these methods will be implemented into design tools.

### (a) High-order mesh generation

Many production simulations with a second-order method require tens or hundreds of millions of cells to produce results of engineering accuracy [6,58]. Some mistakenly believe that high-order simulations would need meshes of similar size. Because high-order methods took much more CPU time than low-order methods on the same mesh, high-order methods were dismissed as prohibitively expensive. In reality, high-order methods can achieve similar accuracy on a much coarser mesh than low-order methods. Therefore, meshes with only tens or hundreds of thousands of elements may be adequate for a high-order simulation. For such a coarse mesh, it is critical to represent curved boundaries with high-order polynomials [33] to achieve high overall accuracy in the simulation.

High-order mesh generation poses two new challenges. First, it is more difficult to generate coarse meshes for a complex geometry, as automated mesh generation algorithms can break down when generating surface meshes at regions with high curvature. Second, generating highly clustered viscous meshes near a curved wall is daunting, as interior mesh lines can cross the curved boundary, or intersect each other, as shown in figure 5*a*. Curved interior elements are necessary to remove the crossing, as shown in figure 5*b* for quadratic triangular elements.

Several approaches have been used to overcome some of the difficulties. Many groups generated fine multi-block structured meshes first. Then these fine meshes are merged once or twice to produce high-order quadratic (degree 2) or quartic (degree 4) quadrilateral and hexahedral meshes. An example of merging four quadrilateral elements into one quadratic element is displayed in figure 6. Although this approach enabled high-order simulations to be carried out, it is time consuming to generate structured meshes for complex geometries and is not, therefore, a long-term solution. Another approach is to generate a linear mesh as coarse as possible using commercial mesh generators. Then the elements with a curved wall boundary are made high order by generating curved edges and surfaces, as shown in figure 5*a*. Finally, the interior elements are also curved based on solid mechanics [59] to avoid grid lines crossing into each other, thus guaranteeing the positivity of the Jacobian of the geometric transformation. In case the geometry is not available, a surface reconstruction technique is used to rebuild a high-order surface [38,60] before the surface and volume meshes are curved.

Ultimately, aircraft manufacturers as large as Boeing or Airbus may not be able to develop and maintain its own high-order mesh generator. Commercial high-order mesh generators need to be developed for these high-order methods to be used routinely in design.

In order to facilitate communications between a high-order mesh generator and a high-order flow solver, it is necessary to define a standard format to store a high-order hybrid mesh containing the usual types of elements, including triangular, quadrilateral, tetrahedral, hexahedral, pyramidal and prismatic elements. A couple of years ago, a public domain mesh generator named Gmsh [61] was the only one capable of supporting high-order elements, and Gmsh format was therefore selected by the First International Workshop on High-Order CFD Methods. Recently the widely used CFD standard called CGNS [62] was successfully extended to handle cubic elements, and a new proposal supporting quartic elements was approved. The existence of such a mesh standard for high-order elements is very critical, as there are so many more variations in high-order elements than linear elements.

### (b) Shock capturing

Shock waves have to be dealt with in the design of future high-speed transonic and supersonic aircraft. Shock-capturing low-order finite volume methods have demonstrated their capability in aircraft design. For adaptive high-order methods, there are two main approaches: limiter [34,63,64] and artificial viscosity [65,66]. There are pros and cons to each approach, and neither is fully satisfactory. The ultimate shock-capturing approach should satisfy all the following requirements:

— accuracy preserving away from the discontinuity;

— free of user adjustable parameters;

— capable of converging to machine zero for steady problems; and

— positivity preserving for pressure and density [67].

Although currently no approach satisfies all the requirements, there has been sufficient progress in both approaches to allow high-quality steady and unsteady simulations with shock waves to be performed. It appears that the limiter approach can be made essentially parameter free, and accuracy preserving, but is often difficult to achieve iterative convergence for steady problems, whereas the artificial viscosity approach is convergent, but not parameter free. Generally speaking, limiter is preferred for unsteady problems, and artificial viscosity is more popular for steady problems. For example, a positivity and accuracy preserving compact WENO limiter was developed and applied successfully to various unsteady problems with strong shock waves in a very robust manner [67], including the double Mach reflection problem shown in figure 7. Several artificial viscosity-based approaches were successfully employed to compute steady inviscid, viscous and turbulent flows with shock waves [68,69].

Shock-capturing methods degrade to first-order accuracy locally near a discontinuity because the error in the location of the shock is proportional to the mesh size. Methods that offer natural subcell resolution can make the error smaller, but cannot change the order. This argument suggests *h*-refinement near shock waves, coupled with a piecewise constant reconstruction, which is the robust, first-order Godunov method. How a locally first-order scheme affects the solution elsewhere is not clear, especially for unsteady flow problems. If the mesh near a shock wave is sufficiently fine such that the magnitude of the first-order error is comparable to the high-order error elsewhere with a coarser mesh, the first-order method will clearly not affect the overall accuracy of the simulation. In any case, *hp*-adaptations for problems with shock waves appear to offer the best promise in accuracy and robustness.

### (c) Efficient time integration algorithms

When the DG method was first developed, the Runge–Kutta algorithm was used for time integration [35]. For unsteady flow problems on a relatively uniform mesh, it was quite adequate. For high-Reynolds-number turbulent flow problems, the resolution of the viscous boundary layer demands an anisotropic mesh near solid walls with aspect ratio as large as 10 000. These anisotropic elements impose an extremely severe time-step limit for any explicit schemes and make them ineffective. Therefore, implicit time marching schemes are essential for such problems. For RANS turbulent flow simulations, the highly nonlinear nature of turbulence models introduces extra stiffness into the already very stiff high-order operators, making convergence to the steady state very difficult and time consuming.

There has been significant progress in the development of time integrators for high-order spatial operators in the past decades, e.g. the LU-SGS algorithm [70,71], *hp*-multigrid solvers [72–75,76], Krylov subspace methods such as GMRES with various preconditioners [77,78] and mixed explicit/implicit approaches [79]. The most used turbulent flow solvers appear to employ a GMRES algorithm with either an ILU, *hp*-multigrid or a line preconditioner [75]. Techniques to improve the robustness of turbulent flow computations such as line search were demonstrated in [80]. With this progress, it is now possible to obtain three-dimensional steady RANS turbulent flow solutions with little user interference.

Shown in figure 8 is such a turbulent flow computation at a transonic speed over a delta wing with artificial viscosity for shock capturing and the *k*–*ω* turbulence model [69]. Key flow features such as the shock wave and vortices were resolved well with a fourth-order DG method.

For RANS simulations, high-order methods are still not robust and efficient enough to be used as a design tool. The stiffness can still drive simulations divergent, or stall convergence to a steady state. In addition, the memory requirement grows with the order to the power of 6 in three dimensions. This is because *n*_{df} is proportional to *p*^{3}, and then the matrix size is of the order of *n*^{2}_{df}. Memory may become a bottleneck for very high-order space discretizations, i.e. *p*>3. Fortunately, for most real-world problems, it appears *p*=2 and 3 offer the most benefit as shown in figure 4. Research is very much needed to improve the robustness and the rate of convergence and reduce the memory requirement of implicit solvers. Finally, the scalability and performance of these implicit approaches on massively parallel computers is another topic of considerable interest.

### (d) Error estimate and *hp*-adaptations

Mesh adaptation (or *h*-adaptation) has been demonstrated in academia and government research laboratories for low-order methods as a very effective way to reduce simulation cost and improve solution accuracy, especially for unstructured grid-based tools. Its use in commercial CFD tools has not been as widespread as expected probably because of the difficulty in coupling geometric modelling, mesh generation, error estimate, mesh adaptation and flow simulation. Once this software engineering problem is solved, the CFD market will embrace a robust *h*-adaptation tool with open arms.

With the development of high-order methods, the order of accuracy is not fixed to first or second order any more. A typical high-order solver can easily incorporate first- to fourth- or even sixth-order schemes in the same simulation. In theory, one would use a first-order scheme near a shock on a fine mesh, and a fourth-order scheme in a smooth flow region on a coarse mesh. Order adaptation (or *p*-adaptation) adds another dimension into the simulation process. Performing *hp*-adaptation will enable potentially much higher pay-off than having *h*- or *p*-adaptation alone. The decision on where to perform *h*- or *p*-adaptation is not an easy one. In many cases, it is difficult to distinguish a smooth feature from a discontinuity on a mesh with finite resolution. Fortunately, if one does *h*-adaptation instead of *p*-adaptation with a smooth feature, no serious harm is done to the simulation.

There has been great progress in the past decades on error estimates and *hp*-adaptations, and interested readers can refer to a recent review [81] for more details. Here, we highlight the success and potential demonstrated by the goal-oriented adjoint-based error estimate and adaptation [82–86]. Let *Q*_{h} denote the numerical solution, and *J*_{h}(*Q*_{h}) the scalar engineering output of interest (such as the lift or drag coefficient). The output adjoint *ψ*_{h} has the same dimension as *Q*_{h} and is the sensitivity to an infinitesimal residual perturbation, which can be computed based on the discrete adjoint equation
4.1
The adjoint has become a very powerful tool in error estimation and adaptations. Its application to high-order methods appears to speed up CFD simulations significantly, sometimes by orders of magnitude.

Next, we present a sample application of adjoint-based *h*-adaptation [87] for turbulent flow over a flat plate at *M*=0.2 and *Re*_{L}=10×10^{6}, a benchmark case from the Turbulence Modelling Resource website. The drag on the flat plate was used as the goal. The DG computations were performed with MIT's ProjectX (PX) code. Computational results from well-known NASA second-order finite volume codes, CFL3D (structured grid-based) and FUN3D (unstructured grid-based), are included for comparison. In addition, non-adaptive DG simulations are also included. Figure 9 shows *C*_{D} versus . We can observe the following from figure 9:

— Even on the structured meshes provides on the website, the DG discretization gives less error for the same

*n*_{df}than the second-order finite volume results.— The adaptive results are significantly superior to the structured mesh results, and in particular the structured mesh finite volume results. Starting at around 5000 d.f. (corresponding to about

*h*∼0.014), the adaptive results are more accurate than any of the structured mesh results even those of the highest d.f. finite volume cases (with about 160 000 d.f.).

After that, we show another example demonstrating the adjoint-based approach for an unsteady moving boundary problem: two aerofoils pitching and plunging in series at low Reynolds number [83]. The simulation started at *t*=0 and concluded at *t*=7.5. The output of interest is the lift on the second aerofoil integrated from time *t*=7.25 to *t*=7.5. The spatial meshes from the final output-based *p*-adaptation are shown in figure 10. Note that the near-aerofoil and vortex shedding regions are targeted for adaptation, as well as the group of large elements surrounding the mesh motion regions. Comparing with uniform *p* results, the *p*-adaptation was able to reduce the number of d.f. by several orders of magnitude.

It is expected that the relative savings in three dimensions with *hp*-adaptations will be greater than in two dimensions. Impressive results have already appeared in the literature, e.g. in [88]. However, the challenges of coupling error estimates, *hp*-adaptation, mesh generation, geometry and flow solver will be greater too. In addition, the robustness of three-dimensional high-order viscous mesh generation for complex geometries may be the main bottleneck. All these challenges remain topics of considerable research.

## 5. Conclusion

We first review the design criteria for future aircraft. In addition to performance and affordability, environmental impact will become increasingly important. Based on forecasted future growth in aviation, reducing fuel burn, GHG emission and noise become imperative. The US government has established aggressive goals in aircraft performance, fuel burn, GHG emission and noise, which is an important first step. Realizing these goals hinges on new ideas in aircraft and engine technology, and technical breakthroughs in the coming decades. One enabling breakthrough will be high-fidelity simulation tools for aircraft aerodynamics, engine and noise computation.

We believe new generations of design tools for aircraft and engines will be based on adaptive high-order methods capable of handling complex configurations. Preliminary two- and three-dimensional computations documented in the first two International Workshops on High-Order CFD Methods demonstrated the potential of these methods for orders of magnitude improvement in accuracy/efficiency over existing lower-order methods. After a decade of very intensive research in the USA, Europe and other parts of the world, we are inching closer to a new generation of CFD design tools based on these high-order methods. In order of importance, I believe progress in the following areas will make that vision a reality:

—

*Commercial quality high-order mesh generation tools*. Right now, each research group has its own tool to produce high-order meshes. This is adequate for research purposes, but not efficient for production simulations. Robust high-order mesh generators are needed to push these high-order methods into the design process.—

*Robust error estimates and hp-adaptations*. An automated*hp*-adaptation tool can reduce simulation cost by orders of magnitude. Research in this area may produce the greatest long-term pay-off.—

*Highly scalable, efficient, robust and low-memory implicit solvers*. Massively parallel computers including GPUs will be used in future design computations. The scalability of implicit time integration algorithms including preconditioners will be critical.—

*Parameter free, accuracy preserving and convergent shock capturing*. Future transonic and supersonic aircraft, and high-performance engines, demand robust shock capturing.

## Funding statement

The author's research on high-order methods was/has been funded by AFOSR, NASA, DOE, US Navy, NSF, DARPA, as well as Michigan State University, Iowa State University and the University of Kansas.

## Footnotes

One contribution of 13 to a Theme Issue ‘Aerodynamics, computers and the environment’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.