{'fraction_non_alphanumeric': 0.029908523512315974, 'fraction_numerical': 0.008635475725379961, 'mean_word_length': 4.962431645479724, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Industry-Relevant Implicit Large-Eddy Simulation of a High-Performance Road Car via Spectral/hp Element Methods

We present a successful deployment of high-fidelity Large-Eddy Simulation (LES) technologies based on spectral/hp element methods to industrial flow problems, which are characterized by high Reynolds numbers and complex geometries. In particular, we describe the numerical methods, software development and steps that were required to perform the implicit LES of a real automotive car, namely the Elemental Rp1 model. To the best of the authors' knowledge, this simulation represents the first fifth-order accurate transient LES of an entire real car geometry. Moreover, this constitutes a key milestone towards considerably expanding the computational design envelope currently allowed in industry, where steady-state modelling remains the standard. To this end, a number of novel developments had to be made in order to overcome obstacles in mesh generation and solver technology to achieve this simulation, which we detail in this paper. The main objective is to present to the industrial and applied mathematics community, a viable pathway to translate academic developments into industrial tools, that can substantially advance the analysis and design capabilities of high-end engineering stakeholders. The novel developments and results were achieved using the academic-driven open-source framework Nektar++.

1. Introduction. The accurate modelling and prediction of fluid dynamics structures and aerodynamic forces around complex bodies is a classical and challenging problem that lies at the heart of modern automotive design. The increasing availability and lowering cost of computational modelling compared to experimental testing means that computational fluid dynamics (CFD) now forms a key part of the engineering design process. When considering the high Reynolds number flows around the complex bodies found in these problems, there are two key considerations that need to be taken into account: (a) creating highly accurate meshes for the geometry at hand, and (b) producing physically-faithful results, while using numerical discretisations that are necessarily unable to capture all the fluid dynamics scales involved due to the large computational costs required.
In industry, the most common practice revolved around the use of ensembleaveraged steady-state modelling, commonly using either Reynolds-averaged Navier-Stokes (RANS) or, more recently, Detached eddy Simulation (DES). Typically these are used in conjunction with low-order numerical discretisations, predominantly the second-order accurate finite volume method, together with linear meshes of planarfaced elements. Although this approach is computationally cheap, it can undermine the accuracy of the simulations and keep the practitioner from identifying aerodynamic design solutions that would be otherwise feasible if the CFD methodology was able to more faithfully describe the flow physics [75]. Indeed, on one hand, the use of highly diffusive numerical schemes and steady-state modeling strategies leads to suppress under-resolved scales and their feedback on the mean flow, thereby additionally corrupting the accuracy of the simulation. On the other, the use of linear meshes requires high degrees of refinement in order to capture highly-curved surfaces in the geometry, which may introduce artificial diffusion and significantly alter the flow evolution that one would otherwise observe. Addressing point (a) and (b) within the context of industrial applications is therefore of primary importance to advance the state-of-the-art of CFD, and can dramatically enhance the computational aerodynamic analysis and design capabilities of high-end engineering stakeholders.
Point (a) requires substantial theoretical and algorithmic advances in high-order mesh generation. From this perspective, the aim is to produce a methodology for generating high-order meshes that achieve levels of robustness over a range of very complex cases comparable to that found in commercial linear mesh generators. Indeed, commercial packages, such as Pointwise [66] and Centaur [8], have attempted to add high-order meshing capabilities to their products with variable degree of success.
The standard approach in generating a high-order mesh is to adopt an a posteriori procedure, whereby a coarse linear mesh is first generated, and is then deformed by introducing the geometric curvature into elements lying on the surface. Several research groups have been working on such an approach, see for example [13,74,44,72,70,85], which gives an indication towards the level of interest in this area.
The challenge of this approach, however, lies in the second step of the process: introducing curvature purely at the surface will typically lead to elements becoming inverted and unsuitable for computation. This has proven a significant stumbling block in the use of high-order methods in industry, and many of the aforementioned works focus on developing techniques to regularize or heal the mesh as part of the a posteriori process. While this is the main route to high-order mesh generation, some alternatives exist. These include the use of the Nash theorem to map a Riemannian space onto an isometric Euclidean space of higher dimension [62,87,11], semi-structured approaches to block mesh generation [3], anisotropic high-order boundary-layer mesh generation [61], and PDE-based approaches based on elastic analogy [68].
In this work, we tackle the problem of high-order mesh generation by using a variational framework that uses the minimization of an energy functional. The process begins by defining a linear surface mesh that conforms with the underlying geometry, typically defined in a computer-aided design (CAD) format that follows standard ISO 10303 (informally known as standard for the exchange of product model data, or STEP) [69]. The linear surface mesh is then complemented with a linear volume mesh along with a boundary-layer mesh (if required). Finally, the linear surface mesh is made high-order to fit the curved surfaces of the underlying geometry, and the interior mesh entities (edges and surface) of the linear volume mesh are deformed to accommodate the curved boundaries imposed by the geometry. The latter step is achieved by the minimization of the energy functional mentioned above, that can be arbitrarily selected by the practitioner depending on the mesh requirements. The linear mesh is obtained beforehand, and uses a combination of the commercial packages CADfix and Star-CCM+ for CAD healing and for linear and boundary-layer mesh generation, respectively, as well as NekMesh [80], which is the open-source high-order meshing framework developed as part of the Nektar++ project [63,5,59].
The second point, (b), is instead closely linked to the impossibility of carrying out simulations that can resolve all the scales present in the problem, that, in turn, leaves part of the flow physics marginally or under-resolved. The computational power required to perform a fully resolved flow simulation, also known as direct numerical simulation (DNS), is a function of the Reynolds number [37]. From this perspective, we can consider the number of floating point operations per second (FLOPS) that a computer can perform as an approximate measure of computational power, and we can then model the time to solution under different regimes of computational power growths. This gives a rough estimate of the number of years required to achieve DNS capabilities for different Reynolds numbers. Figure 1 depicts two different computational power growths. The red line assumes doubling flops every three years, while the blue line assumes a doubling of FLOPS every one and a half years. They represent a less and more optimistic scenario than Moore's law, respectively. We can appreciate how, under the most optimistic power growth scenario (blue line), DNS capabilities for automotive simulations, which corresponds to Reynolds numbers of the order of 10 6 , will be reached approximately by 2040. Figure 1 is an ideal scenario, as it does not take into account additional costs incurred, and it also assumes using the top supercomputers, as soon as they become available. The latter aspect is obviously not achievable for analysis and design purposes in an engineering workflow. Therefore, the ability to perform DNS for industrial purposes at high Reynolds numbers on a regular basis is, at least, two or more decades ahead, even with current advances in computational technologies (e.g. many-core architectures, graphical processing units, and quantum computing). Given these limitations, we need to deal with under-resolution for the next several years, while pushing the boundaries of high-Reynolds number flow physics for industrial problems. Computational resources required to perform the direct numerical simulation (DNS) as a function of the Reynolds number [37]. Yellow and red lines represent projections on when DNS capabilities are due to become available as computational power increases following two different Moore law's paths.
A significant body of work has been produced in the past several years to deal with high-Reynolds number problems. As briefly mentioned above, RANS and DES are frequently used in the context of industrial CFD, see for example [23,86,4], but they typically fail to capture unsteady flow features, which may be crucial for the analysis and design of a given aerodynamic solution. This aspect is especially marked when these two modeling choices are adopted in conjunction with low-order numerical schemes, that are commonly characterized by significant diffusion and dispersion, thereby causing an unphysical behaviour of the flow physics [7].
In order to overcome these issues, several research groups have been working on LES, both physical (or explicit) and implicit, over the past few decades. Explicit LES uses physical principles such as eddy viscosity to devise closure models that can approximate the under-resolved scales, while implicit LES utilizes the diffusion introduced by the numerical discretization. For a detailed review on LES methods, the interested reader can refer to Pope [67]. LES is particularly promising for automotive applications, since these regimes heavily feature flow physics that steady-state approaches struggle to capture, such as vortex separation and interaction [75]. In addition, when used in conjunction with low-diffusive and dispersive numerical schemes, the resolution power is enhanced and can produce flow simulations at unprecedented accuracy.
In this work, we use an implicit LES approach, where we treat the under-resolved scales via a novel artificial numerical regularization approach based on a so-called spectral vanishing viscosity (SVV) kernel, and we adopt robust dealiasing techniques to mitigate aliasing driven numerical instabilities that are normally present due to the low-diffusive behavior of the numerical framework. In particular, we show that the use of SVV and dealiasing strategies, along with key improvements in the computational efficiency of kernel calculations, leads to a robust (in terms of numerical stability) and accurate high-fidelity LES framework for the analysis and design of complex industrial problems. This is showcased in the successful high-fidelity simulation of the real geometry of the Rp1 automotive car at a Reynolds number of 10 6 based on the car length. All the simulations performed use the incompressible Navier-Stokes solver implemented in the open-source spectral/hp element framework Nektar++. The solver implements an operator splitting scheme [36,29], where the spatial discretization is achieved via a continuous Galerkin method [38], and the timediscretization is obtained via a second-order implicit-explicit time-stepping strategy [83]. The numerical simulations performed were achieved with hybrid meshes containing both prisms (to discretize the boundary-layer regions in proximity to the body of the car) and tethraedra (to discretize the farfield of the domain).
In this paper, we present the key advancement in high-order meshing and simulation technologies that we implemented in NekMesh and Nektar++ in order to achieve the challenging simulation of the Rp1 car, with the aim translating academic research efforts to industrial CFD applications. In particular, section 2 describes the high-order meshing strategy we adopted, focusing on the steps required for surface and linear mesh generation as well as on the a posteriori high-order deformation. In section 3, we describe the novel simulation technologies implemented in terms of regularisation approaches for the under-resolved scales, dealiasing techniques for numerical stability, and computational efficiency. Section 4 presents the high-fidelity simulation of the Rp1 car, and finally, in section 5, we draw some final remarks.
2. High-order mesh generation for complex geometries. The starting point and essential prerequisite for any high-fidelity LES, within either a finite el-ement or finite volume framework, is the generation of a geometrically-accurate and high-quality mesh. Mesh generation is a complex and involved process, particularly when the geometry under consideration is a complex, three-dimensional body, since the important regions that need to be resolved within the domain are often unknown a priori. The process may therefore require many iterations of meshing, simulation and analysis in order to arrive at a final mesh deemed suitable enough for accurate results.
Mesh generation tools are generally required to be both flexible and robust, leading to a wide variety of approaches. In industry, the majority of mesh generation strategies lead to a linear mesh, as this is what most simulation packages are capable of supporting. While widely adopted, the use of linear meshes could alter the flow physics in proximity of curved geometrical features, that in turn can lead to a critical reduction in simulation accuracy. The use of higher-order methods can address this loss of accuracy, as one generates boundary-conforming, curvilinear meshes which describe more accurately the underlying geometry for the purposes of simulation.
However, the generation of such meshes in a robust and efficient manner for a realistic three-dimensional geometry remains an open and challenging problem. In fact, the introduction of additional curvature within elements that are connected to the geometry -and that must conform to the boundary -will often lead to significant reduction in the quality of the mesh, or indeed deform the element to such an extent that it is self-intersecting and thus unsuitable for simulation usage. For simulations in fluid dynamics, where boundary layers of highly anisotropic elements are required to resolve the shear of the flow near the boundary, this issue becomes even more pronounced. Finding routes to mitigate and address these issues is the central challenge for the high-order mesh generation community, and forms a key bottleneck to tackle in the wider adoption of high-order methods. In this section, we outline our process to try to overcome this issue, and highlight the process within the context of realistic automotive geometries.
2.1. Approaches to high-order a posteriori mesh generation. Our framework to produce a viable high-order mesh, similar to other state-of-the-art methods [13,74,70,32,85], uses an a posteriori approach. In this setting, one first generates a coarse, linear mesh of the desired geometry. From this, high-order nodes are added within each element. Any nodes that lie on curved surfaces are projected onto the surface in such a manner as to enforce a boundary-conformal representation. In order to achieve the curvature required by the geometry being meshed, we typically adopt an isoparametric approach, in which the same basis functions that are used to represent high-order solutions within each element are also used to define the geometric coordinates within the element. Such a mapping can be expressed therefore as the summation where x is the spatial location of the a point within an element, ξ is a corresonding location within a reference element, x n is the location of a given node in the element, and n is the polynomial interpolant (e.g. a Lagrange polynomial interpolant). In particular, mesh entities, namely faces and edges, on the boundary of the geometry are curved by adding high-order nodes, and then enforcing that the nodes lie on the boundary. Moreover, the position of the nodes x n must be adjusted so that the resulting mapping is a faithful polynomial representation of the geometry [74].
Although from this description this may seem a simple task, for complex geometries that can have several hundred curved surfaces, the ability to create valid high-order elements (that is elements of order greater than 2) is extremely challenging. In many cases, when simply curving the edge and/or face in proximity to a given curved geometrical feature, it is relatively easy to generate invalid elements. Figure 2 shows seven mesh elements in proximity of a curved surface, depicted in grey. The first mesh element (a) is a valid linear element that completely misses the curvature of the underlying surface. Indeed, in order to capture the curvature of the surface with a linear mesh it is necessary to add more elements around the surface than would be necessary with high-order elements; the latter benefits from a superior accuracy that allows for larger element sizes than linear elements. The second element (b) is a valid high-order element, where only the edge in proximity to the geometry is curved.   that the challenge of high-order mesh generation is multifold. On the one hand, it is necessary to create curved mesh entities that conform to the underlying curved geometrical features (where the geometry is typically provided in a CAD format).
On the other hand, it is essential to develop a strategy to handle invalid high-order elements and produce a valid high-quality mesh. Our high-order meshing strategy handles the challenges just described and consists of three steps, (1) the definition of the geometry in such a way that is can be efficiently queried for curvilinear mesh generation, (2) the generation of a highquality linear mesh, and (3) the generation of curved elements. These three crucial steps are described in detail in sections 2.2, 2.3 and 2.4, respectively. The corresponding meshing framework has been implemented in the code NekMesh, that is part of the open-source Nektar++ library [5]. For a more comprehensive explanation of the mesh generation process, the interested reader can refer to Turner et al. [81].
2.2. Geometry definition. The first stage of our meshing approach is the definition of the geometry in a way that ensures its accurate representation and permits geometrical enquiries during the various mesh-generation stages. The most common approach to define an arbitrary complex geometry is by a boundary representation (also referred to as "BRep"). The boundary representation contains surface patches and curves that link these surface patches. The surface patches and curves defining the underlying geometry are typically topologically connected, thereby forming a closed shell. These surfaces and curves are parametrised onto a lower-dimensional space, that is the one-dimensional space x(s) for curves, and the two-dimensional space x(s 1 , s 2 ) for surfaces, where x = (x, y, z) denotes the three-dimensional coordinate of a point in the boundary representation of the geometry. The main geometrical enquiries required for mesh generation consists of three key operations: • parametric coordinate to real space location, i.e. determining x from (s 1 , s 2 ); • projection of a location in space onto a boundary representation entity, i.e. determining (s 1 , s 2 ) from x; and • calculation of the first-and second-order derivatives of the parametric mapping with respect to s 1 and s 2 . In addition, in this work we do not consider the effects of incomplete or low-quality geometric representations; The boundary representation must be watertight and should represent accurately the geometry under investigation. These two requirements are crucial for high-order mesh generation, as they impact the number of high-order valid elements that is possible to generate and the quality of the mesh.
The assumption adopted in this work is that the geometry is given in a standard CAD (computer-aided design) format that supports a boundary representation description, such as STEP. In several CFD applications, this assumption is fairly common, especially among industry practitioners. The CAD object is then probed via NekMesh, that performs geometrical queries through interfaces to a CAD engine. To give support for different CAD engines and implementations, NekMesh implements a lightweight interface to the key operations above. Concrete implementations of this use CADfix [34] and OpenCASCADE [64], both of which permit reading CAD data in various formats, including STEP. In particular, the checks and queries on the CAD object were performed by NekMesh, while the correction of any errors was performed by CADfix which provides an extensive toolset for fixing geometrical inconsistencies in CAD data. CADfix is a commercial software, that provides extensive CAD handling tools, useful for healing and modifying the underlying CAD geometry definition. It also provides a linear mesh generator based on the medial approach that is useful for decomposing the domain into partitions. This decomposition has proven useful for designing high-quality boundary-layer meshes near the wall surfaces, aspect that was critical for the high-fidelity simulations in section 4. OpenCASCADE is an opensource code with reduced capabilities than CADfix, and represents the default CAD handler within NekMesh, due to its free-availability, that is at the core foundations of the NekMesh and Nektar++ projects. Nevertheless, the use of CADfix in this work was essential to obtain a robust and efficient high-order meshing pipeline, and ultimately obtaining valid high-quality meshes. Access to CADfix, or some alternative approach to geometry healing and defeaturing, is a fundamental requirement for the development of a robust and successful meshing strategy.
Once the CAD representation of the geometry is satisfactory, it is possible to complete the definition of the computational domain, including the external boundaries, and proceed to the next meshing stage, that is generating a valid linear mesh. This step is described next.
2.3. Linear mesh. The generation of a high-quality linear mesh from the CAD representation of the geometry is a demanding task. Indeed, it is necessary to have a sufficiently coarse initial mesh that serves as a starting point for the high-order meshing process. This initial mesh is typically obtained via a user-driven manipulation of the mesh spacing definition at all points within the computational domain to achieve satisfying geometric accuracy, coarseness and mesh gradation. This procedure is typically time-consuming and challenging, possibly involving a large number of iterations before obtaining a linear mesh that satisfies the quality required by the practitioner.
In order to automate this step as much as possible, we combined the advanced linear meshing capabilities of the commercial package Star-CCM+ with the high-order meshing procedures of NekMesh. This allowed us to produce an automatic industrial linear meshing pipeline, that required minimal intervention from the user. The main obstacle that we encountered when using Star-CCM+, and that is common to encounter when using commercial CFD software, was the lack of information relating the linear mesh to the original CAD definition. Often the CAD surfaces are replaced by a fine triangulation that represents the "true" boundary of the geometry within a given tolerance. The linear mesh is then built using this approximate representation and the original CAD information discarded. This is a significant issue since inaccuracies in the definition of the geometry could lead to errors in the high-order mesh generation process, that in turn could lead to inaccurate flow simulation results.
In order to tackle this problem, we used the projection of the linear mesh onto the CAD surfaces to reconstruct the missing CAD information. In particular, we imported the linear mesh generated by Star-CCM+ into NekMesh. The import was performed by incorporating curvature-driven mesh size restrictions such that the boundary of the surface mesh is a reasonably close representation of the underlying CAD model. Once the import of the linear mesh into NekMesh is complete, the CAD model is processed in order to obtain parametric coordinates in the CAD representation that will be used to enhance the adherence of the linear mesh to the CAD model. From this perspective, the first step is to linearise each surface of the CAD model via an "auxiliary" triangulation, so that each CAD surface has two representations, the CAD model itself, and the "auxiliary triangulation". The second step is to define a bounding box (inflated by 5% in each direction), for each CAD surface, and store these boxes into a k-dimensional (k −d) tree data structure [1]. Once the two CAD representations are available and the bounding boxes are defined for each CAD surface, the surface mesh nodes are processed to obtain their CAD parametric coordinates. More specifically, we associate each node of the surface mesh that was generated in Star-CCM+ to a list of potential CAD parent surfaces by querying the bounding box k − d tree. If a node is within a bounding box, the surface is added to the list of potential parents surfaces. The use the k −d tree makes this process particularly efficient and reduces the number of potential parent surfaces to just a few candidates, typically on the order of 2 or 3. The surface mesh node is then projected onto the CAD surface for each of the parent surface candidates. This is achieved by identifying the nearest vertex in the "auxiliary triangulation" of the CAD surface, and then using this vertex as an initial guess for the nonlinear optimization problem associated to node-to-surface projection. The parent surface is finally identified as the closest surface to the projected node.
We note that, since the "auxiliary triangulation" of the CAD surfaces is not directly related to the Star-CCM+ linear mesh, the surface node of the linear mesh might not lie on the CAD surface. Therefore, once the parent CAD surface associated to the mesh node is identified, the mesh node is moved to the CAD surface it belongs to, thereby improving the adherence of the linear mesh to the underlying CAD surfaces. The node movement just described is always performed, except in two cases: if the required displacement of the node is greater than 10% of the length of the edges in the local mesh, and if the movement of the node induces inverted (i.e. invalid) elements. In these two instances, the node is added to an "exclusion set" with a view to preventing deterioration of mesh quality or generation of invalid elements. The final product of the projection process just described is a set of nodes that are accurately placed onto the original CAD surfaces, and a set of nodes that are placed into an "exclusion set". Each node, regardless if in the "exclusion set" or not, will have all the parametric information required to generate the final high-order curving of the surface mesh, with the exception that any mesh entity (edge or face) that has a node in the "exclusion set" will be left linear. The approach outlined has been named projection curving [80], as it involves the projection of the surface mesh generated via Star-CCM+ onto an "auxiliary triangulation" of the CAD surfaces. This approach is the one used in the final industrial meshing pipeline, that was adopted for the last design stage, D3, of the road car in section 4. A boundary-layer mesh was also created by Star-CCM+, and consisted of a prism layer in proximity of the surfaces of the geometry.
We finally point out that we initially attempted two additional alternatives, which were however discarded because they became too time-consuming and user-driven. These are named analytic curving, and were used for the first two designs, D1 and D2, of the road car presented in section 4. In both the analytic curving strategies, all the parametric information associated with the linear surface mesh was known, hence high-order curving of the surface was a relatively simple task, unlike in the projection curving approach. The first analytic curving strategy used the linear meshing capabilities of NekMesh to create the surface mesh that was later exported to Star-CCM+ for the generation of the volume mesh, and of the boundary layer mesh (constituted by a near-wall macro prism layer). This approach was however extremely time-consuming as it involved dozens of iterations between NekMesh, Star-CCM+, and the original CAD model. The second analytic curving strategy tried to address the bottlenecks of the first by allowing the Star-CCM+ linear mesh generator to create the surface mesh from a linearised CAD surface obtained using NekMesh. The surface mesh was then forced to explicitly conform to the boundaries of the CAD surfaces. This is critical in an industrial meshing pipeline, as it is typically impractical to clean complex CAD models with thousands of surfaces, so as to guarantee a clean enough CAD to enforce the surface mesh to obey to the CAD surfaces.
2.4. Generation of high-order elements. The third stage is the generation of the curved surfaces, where the coarse linear mesh generated in the previous step and described in section 2.3 is deformed to accommodate the curvature of the boundary. We already mentioned in section 2.3, that the creation of curved boundary mesh entities that conform to curved geometric features can be handled via two different strategies, projection curving and analytic curving, and we found the former to be the optimal approach for complex industrial CAD models. While curving the boundaries attached to the geometry is a relatively simple task once the parametric coordinates of the surface mesh are known and accurate, the creation of valid high-order elements remains a challenge, as depicted in figure 2. Indeed, one of the key research themes in a posteriori high-order mesh generation is to find suitable deformations of the interior mesh entities. The ability to deform interior mesh entities allows for better mesh quality and for a significant reduction of invalid elements, see figure 2-(e). To obtain suitable deformation of interior mesh entities, there exist a certain number of alternatives, including methods based on partial differential equations [65,85,60,22,30] and on the minimization of an energy functional [78,24]. In this work we pursued In this example the problem could easily be solved by curving the other two edges of the triangle. But whereas the domain definition, usually from CAD, defines how to curve the boundary mesh entities, no such guidance is available on how to curve the interior entities. Therefore the task of high-order meshing is to obtain a deformation such that the resulting curvilinear mesh is geometrically conforming with the boundary and its interior is curved to accommodate the boundary deformation and produce a valid mesh. However, it is notoriously di cult to obtain a the latter, and developed a variational framework able to encapsulate arbitrary energy functionals in a unified high-order meshing framework that has been implemented in NekMesh.
The variational formulation begins with a coarse linear mesh that was created as described in section 2.3. Following [80], we consider a linear mesh with N el straightsided elements, that is Ω I = N el e=1 Ω e I , where each element is provided with a highorder polynomial basis function based on standard Lagrange polynomials. The mapping between the straight-sided mesh, Ω I , and the curvilinear mesh, Ω, is denoted by φ : Ω I → Ω.
With reference to figure 3, we denote the coordinates of the standard (also referred to as reference) element as ξ ∈ Ω st , those of the curvilinear element as x ∈ Ω e , and those of the ideal element as y ∈ Ω e I . The mappings between the elements in figure 3 are constructed in an isoparametric fashion, such that the nodes ξ n defining the Lagrange basis functions on the standard element map to y n under φ I and to x n under φ M . The example reported in figure 3 is for triangular elements but it can be readily generalized to any element shape, in both two and three dimensions.
Once the nomenclature for the mesh is specified, we can define the energy functional adopted throughout this work, that is where W depends on the deformation gradient tensor ∇φ(y) = ∂φ ∂y and its determinant J = det[∇φ], also referred to as the Jacobian. The functional W can assume various forms depending on the associated energy adopted. We have implemented functionals associated to linear elasticity energy, isotropic hyperelasticity energy, Winslow equation energy, and a distortion-based energy. Details of these functionals can be found in Turner et al. [81] and in Turner [80], and are not reported here for the sake of brevity.
Since the ideal mesh is a union of many elements, the functional can subsequently be defined element by element where we used the composition φ = φ M • φ −1 I . To fully define equation (2.3), we need to specify the expression of φ M and φ I . The first, φ M , assumes the usual form of an isoparametric mapping commonly found in nodal spectral element methods and specified in equation (2.1). The second, φ I , assumes an analytical form that is a combination of linear finite element basis functions and that depends on the element shape. For instance, in the case of a triangle with vertices v 1 , v 2 , and v 3 , the mapping is given by The minimization of the functional is then formally defined as follows where we seek the configuration of nodes of the linear mesh, denoted by N = [n 1 , . . . , n N nodes ], that minimizes the energy functional in equation (2.5). In practice, this means to start from the initial distribution of nodes in the linear mesh, N 0 , in which the curvature of the boundary mesh entities is imposed, regardless of whether it causes self-intersections. Then, the minimization proceeds in iteratively updating N towards lower energy configurations, that is: The minimization is achieved via a gradient-and Hessian-based method, and utilizes a relaxation strategy, whereby a local optimization problem is solved in place of a global one. More specifically, the minimization in equation (2.5) is performed for a subset of elements that are influenced by a change in the position of the node under consideration. This requires the solution of a smaller and less computationally expensive problems, compared to the global problem that would be otherwise necessary to solve. The minimization of equation (2.5) can then be recast into the minimization of a set of energies, each of which belongs to a node and the subset of elements the node influences. The stopping criterion used to find the global minimum based on the set of local energies, is ||N k+1 − N k || ∞ /L < , where we used = 10 −6 , and ||s|| ∞ = max{|s 1 |, . . . , |s n |} is the uniform norm.
This process also requires the computation of the gradient of E(φ), with respect to the movement of node positions. As a first approach, these gradients can be approximated by use of finite difference stencils. However, our testing indicates that this approach is both expensive and sensitive to the size of the stencil. Instead, we make use of an analytic expression for the desired gradients. Ultimately, applications of the chain rule therefore require the same derivatives for φ I and φ M . Although the ideal mapping φ I readily admits an analytic expression, more care is required for φ M , but by exploiting properties of the Lagrangian expansion (2.1), analytic forms of the derivatives can be obtained, which are given in the appendix of [81]. Figure 4 shows the high-order surface mesh of the Rp1 car, resulting from our meshing process. It is possible to appreciate the high-order elements tessellating the CAD surface. We note that the spokes within the wheel shown in figure 4 were not included in the simulations presented in section 4 to avoid the requirement for a moving mesh. Fig. 4. Details of the Rp1 high-order mesh produced using the high-order meshing workflow presented in section 2. We note that the spokes in the wheel were not included in the final simulations to reduce complexity.
As a final remark, we point out that, all those nodes on the boundary surface mesh that belonged to the "exclusion set" are left straight-sided. This leads to a small but not insignificant proportion of the surface mesh entities that will be straight-sided as opposed to curved, and could have significant impact on the geometric accuracy of the mesh and hence the solution. However, in the results that we present in section 4, we observe little noticeable impact, since often these excluded elements were in locations where multiple CAD surfaces met but not typically in places of interest to the flow physics. Obviously this cannot be always guaranteed to be true. A full investigation into this has not been conducted yet, but this work proves that high-order meshes can be produced on industrial models and sensible results obtained with little meshing effort. Indeed the final projection meshing pipeline only required one execution of the linear mesh generation and NekMesh, with no need for repetition with either system or the CAD.
3. Numerical methodology. The model underlying the simulations carried out in this study is constituted by the incompressible Navier-Stokes equations for a Newtonian fluid where u = (u, v, w) is the velocity vector, p denotes the pressure, ν is the kinematic viscosity, and f is a generic forcing term. This set of equations are particularly suited to accurately model the flow around a road car, as compressibility effects are deemed to be unimportant at the speeds reached. For diffusers and other heated parts of the car, the compressible Navier-Stokes equations are better suited. This was however out of the scope of the current investigation, which was mainly focussed on the external aerodynamics of the car.
The system of equations (3.1) and its numerical discretization has been implemented in the open-source spectral/hp element library Nektar++. The equations are solved via an operator splitting scheme (sometimes referred to as a velocity correction scheme) [36,29] in combination with a consistent boundary condition for the pressure Poisson equation. In particular, following [29], we recast problem (3.1) into a weak pressure Poisson problem by taking the inner product over the solution domain Ω of equation (3.1a) with respect to the gradient of the test function, ∇v, i.e.
where N (u) = u · ∇u and we set the forcing term f to zero. Using the identity , we can enforce the divergence to be zero by setting the last term to zero. By integrating the first, second and last term in equation (3.2) by parts we get the weak pressure equation where ∂Ω is the boundary of the domain, n indicates the outward pointing normal to the domain Ω and the superscript (n + 1) denotes the time step that aim to compute. Note that in equation (3.3) we have also used the identity ∇ · (∇ × ∇ × u) = 0. By using a backward approximation of the time derivative, along with a consistent extrapolation for the non-linear terms N * ,n+1 , whereũ n+1 is an intermediate velocity, γ 0 is a constant, andû is an extrapolation of the velocity, we can finally calculate the pressure at step n + 1. We can then solve for the velocity at step n + 1 by using the pressure at n + 1 just calculated. This leads to solving the following Helmholtz problem The velocity correction scheme just outlined (equations (3.3),(3.5)) is spatially discretized using a high-order continuous Galerkin formulation [38]. Here, one typically expresses the solution in terms of global basis functions Φ k (x) as follows where u k are the weights or coefficients of the polynomial expansion, and we introduced the superscript δ to denote the spatial discretization of u. Equation (3.6) can be decomposed into local boundary and interior function, such that equations(3.3) and (3.5) can be written in terms of element-wise contributions (omitted here for the sake of brevity). The elemental polynomial expansions can be written as where the elemental basis functions φ k (x) used are modified Jacobi or Lagrange polynomials defined on a set of Gauss-Lobatto-Legendre quadrature points, k are the elemental coefficients of the expansion, and we introduced the element Ω e of our domain discretization. The global-wise continuity of the solution required by the continuous Galerkin method is then enforced by an assembly procedure. The interested reader can refer to [36,29] as well as to [38] (Chapter 4 and 8). The semi-discrete form of the incompressible Navier-Stokes equations that is obtained using the velocity correction scheme and the continuous Galerkin method as spatial discretization and that was just introduced is advanced in time via a second-order implicit-explicit time-stepping strategy [83].
In order to successfully achieve the challenging high-Reynolds number LES described in section 4, we needed to significantly enhance the numerics implemented in Nektar++, especially in terms of (i) numerical stability, (ii) regularization techniques for under-resolved scales, and (iii) computational efficiency. These three key building blocks are discussed in the following.
3.1. Numerical stability. Numerical stability was a key aspect that we needed to carefully take into account in order to successfully simulate the road car presented in section 4, at the Reynolds numbers required by industry. To achieve this, we adopted polynomial dealiasing via consistent integration [48], combined with spectral vanishing viscosity (SVV) [35,40,50]. The former helps suppressing aliasing-driven instabilities that are due to inexact numerical integration of the equations, while the latter targets all those flow scales that are necessarily under-resolved in high-Reynolds number simulations. While SVV primary target was not the numerical stability of the simulations, it was necessary the use of both, dealiasing and SVV, to achieve stable (i.e. non-crashing) computations. More details on the SVV adopted are reported in section 3.2, where we discuss its use to accurately treat under-resolved scales. In this section, we focus on polynomial dealiasing that was a key component without which it was otherwise not possible to obtain stable simulations.
Indeed, the main issue when dealing with element-wise integrations typical of spectral/hp methods, is the inexact numerical integration of the discrete terms of the equations. In fact, it is well known that there is a minimum number of quadrature points Q min , that is required to produce exact results for integrals of non-rational polynomials of any given degree. This minimum number of quadrature points depends on the type of numerical quadrature adopted and on the set of quadrature points used. In our case, we employed Gaussian quadrature on a set of Gauss-Lobatto-Legendre quadrature points [38]. This type of numerical integration guarantees that the minimum number of quadrature points Q min necessary to exactly integrate any given order P polynomial u(ξ) ∈ P P up to machine precision is Q min = (P + 3)/2, where P is the space of polynomials. It is then possible to calculate the number of quadrature points necessary to exactly integrate different nonlinearities that might arise in the equations, as reported in table 1. We can observe how Q min increases as the degree of the nonlinearity increases.
Q ≥ 2P + 3/2 Table 1 Number of GLL quadrature points for the GLL quadrature to be exact up to machine precision as a function of the polynomial order of the integrand Let us now consider the incompressible Naver-Stokes equations (3.1), and split them into linear and nonlinear terms as follows: where M is the mass matrix that arises from the numerical discretization, L = −∇p δ + ν∇ 2 u δ + f δ contains the linear terms, N = −u δ · ∇u δ contains the convective term, and superscript δ denotes the spatial discretization. By noting that the convective term is a quadratic nonlinearity, we can identify the Q min required for exact integration. In particular, given that the spatial discretization we employ is a continuous Galerkin method, we are interested in the integration of L 2 inner products of two polynomials (u p , u q ). On the one hand, to calculate the linear terms, it is necessary to calculate [u(ξ)] 2 ∈ P 2P (i.e. the inner product of the solution with respect to the polynomial basis). Therefore, for exact integration, we need Q min = P + 2 Gauss-Lobatto-Legendre quadrature points in one-dimension and a tensor product of points in higher dimensions. On the other hand, to calculate the nonlinear terms, we wish to calculate [u(ξ)] 3 ∈ P 3P , that corresponds to Q min = 3 2 (P + 1) Gauss-Lobatto-Legendre quadrature points. We note that using other sets of quadrature points, Q min changes. For instance, the use of Gauss-Legendre points guarantees exact integration up to Q min = P + 1/2, for any given order P polynomial u(ξ) ∈ P P . While it is typically possible to exactly integrate the nonlinearities N present in the equations when dealing with mesh elements that are straight-sided parallelograms, it is more complicated if not impossible when dealing with curved mesh elements. In fact, in the case of straight side parallelograms, the Jacobian of the mapping that defines the coordinates of the element is affine with constant determinant. For curved elements however, this mapping is an isoparametric polynomial expansion, which when incorporated into integrands leads to an additional source of aliasing error, called geometrical aliasing [48]. Geometrical aliasing enters into the discrete formulation via the inverse of the mass matrix M −1 , as well as from both the linear and nonlinear terms in equation (3.8), and it is typically in the form of a rational polynomial. It is well known that rational polynomials cannot be exactly integrated within the functional space used to approximate the problem. Therefore, numerical stability cannot be formally guaranteed. We also note that geometrical dealiasing could couple with the underlying equation nonlinearities, thereby exacerbating numerical instabilities. While, for problems that employ unstructured meshes with curvilinear elements like the ones described in this paper, numerical stability cannot be guaranteed, it can be substantially improved by adopting dealiasing strategies that mitigate aliasing-driven instabilities.
The dealiasing approach adopted in this work is the so-called Global dealiasing [48], that employs a larger and identical number of quadrature points on each term of the numerical discretization (3.8) of (3.1), that is The subscript Q in each term of equation (3.9) denotes a higher number of quadrature points than required for exact integration of quadratic linear terms, that is Q > P + 3/2. This dealiasing strategy mitigates aliasing-driven instabilities induced by equation nonlinearities, curved mesh elements, and the coupling between equation nonlinearities and curved mesh elements. Exact integration is however not guaranteed, except for an asymptotically infinite number of quadrature points. Hence, marginal aliasing errors can still be present and contribute to numerical instabilities. In this work, we always used a quadrature order consistent with the nonlinearities of the governing equations, often with a safety margin to account for curved elements, that is Q > 3P/2 + 3/2. We note that for simulations that are adequately resolved (e.g. laminar flows), aliasing usually does not affect numerical robustness. However, for underresolved turbulence computations, this aliasing effect leads to a significant build-up of error in high-frequency solution modes and usually causes the simulation to abruptly diverge.
There are additional tools to suppress aliasing-driven instabilities, that have not been explored in this work. These include the use of skew-symmetric forms of the convective term of the underlying equations [2,25,84], polynomial filtering [27,19,21,31], variational multiscale schemes [33,46], and exponential based modal filtering [26]. As mentioned above, we note that SVV could, in principle, also be used to suppress aliasing-driven instabilities [35]. However, the diffusion required would be significantly high, and could suppress important flow features, thereby undermining the accuracy of the simulation. Additionally, because the SVV operator is relatively isotropic, whereas dealiasing is anisotropic, the two do not completely overlap. Hence, a reduced amount of regularisation can be achieved with dealiasing, thereby better preserving the solution accuracy.
The use of the dealiasing techniques described in this section removes numerical integration errors that could lead to numerical instabilities in under-resolved simulations; however, it does not address the under-resolved scales that are ubiquitous in high-Reynolds number simulations. These under-resolved scales, similarly to aliasing errors, can be one of the main drivers of numerical instabilities and solution inaccuracies, and they exacerbate aliasing errors leading to the possible failure of the simulation. The handling of under-resolved scales to prevent numerical instabilities and maintain solution accuracy is described next.

Regularization for under-resolved scales.
In order to treat underresolved scales commonly present in high-Reynolds number simulations, such as the one in section 4, we use a regularization technique based on a novel SVV operator. This new operator has been carefully designed to incorporate insights obtained from recent dispersion/diffusion eigenanalysis [57,49,47,50], and results in well-balanced diffusion characteristics that allow for the accurate simulation of high Reynolds number flow problems. Indeed, when stabilising a high-order formulation, such as the continuous Galerkin scheme adopted in this study, it is critical to achieve an adequate balance between numerical diffusion and solution accuracy, especially for turbulent flow computations.
The general idea behind SVV is that its diffusion can be designed to affect mainly the element-wise polynomial modes of highest order within the spectral/hp element discretization. These polynomial modes comprise the most oscillatory part of the solution and represent the smallest turbulent scales captured by the degrees of freedom employed. The SVV technique was first proposed by Tadmor [76] as a stabilisation strategy for classical spectral methods in Fourier space. Essentially, if the solution in Fourier space is then SVV ammounts to an artificial diffusion term added to the governing equations in the form where ν svv is a constant representing SVV's overall diffusion magnitude, κ is the wavenumber,û κ are the Fourier coefficients of the solution,Q κ are the Fourier coefficients of the SVV kernel that define how diffusion is to be distributed across Fourier modes, and i = √ −1. We note that by settingQ κ ≡ 1, the standard diffusion (Laplacian) operator is formally recovered. Traditional SVV kernel functions grow monotonically from zero, achievingQ κ = 1 only at the highest order mode. For that solution component, the physical kinematic viscosity ν effectively becomes ν + ν svv .
Since Tadmor's original work on spectral methods, SVV has been adapted to continuous spectral/hp element methods and used in particular to stabilise underresolved turbulence computations based on the continuous Galerkin method [35, 40, 41, 43, 71, ?]. In this approach the filtering properties of the SVV operator applied to the continuous Galerkin (also referred to as CG hereafter) method are used inlieu of an explicit turbulence model. The use of techniques originally designed to promote numerical stability in place of a turbulence model is already a common practice, for example, in model-free turbulence simulations via discontinous Galerkin (or DG, hereafter) methods, where upwind-type dissipation is introduced through the Riemann fluxes employed in compressible formulations [12,82,54]. This setting, where artificial numerical features are used in place of an explicit turbulence model, is commonly referred to as implicit LES, and it has been investigated in the past few decades with various degrees of success [45,56,9,10,52].
The key potential issues when modeling under-resolved scales using the diffusion properties of the underlying numerical discretization or through an SVV-like operator are lack of robustness and/or accuracy of the resulting implicit LES framework, and unclear connections between the artificial numerical filters and the flow physics. In this work, we attempted to address these challenges by encapsulating novel results from the so-called eigensolution analysis (ESA) mentioned previously [57,49,47,50] to build a new SVV kernel. In this context, we sought for solutions in which the element-wise coefficients arising from the numerical discretization are related to wave-like solutions of the form through projection (see also [51] (Chapter 2). Initially, we used numerical dispersion and diffusion estimates for spectral/hp element methods in wavenumber/frequency space provided by the so-called temporal ESA [57] to construct an SVV kernel Q that enforces SVV diffusion to be distributed monotonically across wavenumber space for the CG method. We note in fact that monotonic behavior in wavenumber space is not always guaranteed due the lack of one-to-one correspondence between Fourier and polynomial modes 1 . This resulted in what we called "power-law" CG-SVV kernel where p is the polynomial mode index, P is the overall polynomial order, and P svv is a user-defined parameter that, when increased, further confines CG-SVV damping to the smaller scales. This kernel shape is motivated by hyper-viscosity operators [cite the paper that the reviewer mentions by Avrin and Chang]. We note that the operator Q(p) in equation (3.13) is formally different from the kernel defined in equation (3.11), as this arises from the projection of the element-wise coefficients in equation (3.7) onto the wave-like solution in equation (3.12). This kernel also guarantees that raising P increases the resolution power per degree of freedom which is something expected from high-order methods, as discussed in detail in [57]. We highlight that this first study was based on the temporal ESA framework, which is relevant to temporally evolving simulations instead of spatially evolving ones, where inflow and outflow boundary conditions are typically used. While the "power-law" kernel defined in equation (3.13) solved a number of issues, still some potentially undesirable characteristics were present. These undesirable aspects were two-fold. On the one hand, non-smooth behavior of the diffusion and dispersion curves was observed for either too small or too large values of ν svv . This could lead to unphysical representation of the flow physics. On the other hand, the use of the temporal ESA framework did not include spatially evolving flow problems (i.e. problem with inflow/outflow boundary conditions) that are commonly found in external aerodynamics. This could trigger spurious reflections of certain waves when they hit regions of mesh expansion (e.g. along turbulent wakes), eventually leading to numerical instabilities and divergence of the simulation [49,47]. In order to address the first aspect, we focused on the relevant non-dimensional parameter to quantify SVV's relative intensity is the so-called Péclet number, namely Pe = vh/ν svv , that is based on measures of local velocity v and mesh spacing h. The adopted strategy was to maintain the Péclet number fixed by making ν svv ∝ vh. This led to a CG-SVV scaling that is dimensionally equivalent to the upwind-based dissipation of discontinuous spectral/hp element methods, as explained in [57].
To address the second aspect, we pursued a complementary study [50] focused on the so-called spatial ESA, which assumes inflow-/outflow-type boundary conditions and is therefore more relevant to spatially developing problems, as the one encountered in this paper. Spatial ESA had initially been conducted for discontinuous spectral/hp element methods [49,47] and proved insightful in explaining how spurious reflections occurred in regions of mesh expansion depending on the Riemann flux employed. In this framework, theory prescribes that an incoming wave of frequency ω should be associated with a single (and real-valued) wavenumber κ = ω/a, a being the advection velocity. However, in a numerical setting, two complex-valued wavenumbers, κ p (physical) and κ u (unphysical), appear in the analysis. Regarding their imaginary part, when numerical diffusion is present, one observes (κ p ) > 0, which indicates damping of the physical mode during forward propagation, since the solution can be written as (3.14) u where and correspond to the imaginary and real part, respectively. For the unphysical reflected mode, it is instead necessary that (κ u ) be negative, indicating damping during backward propagation (∆x < 0).
When applied to the CG method in [50], spatial ESA revealed that unless special care is taken during SVV design, unphysical reflected waves can have negligibly small (although negative) values of I(κ u ). This allows for spurious waves to survive for a long time while interacting with incoming turbulence, eventually causing numerical instability. Hence, a new optimisation procedure was performed for the CG-SVV operator in order to match the diffusion levels of CG to those of standard upwind DG (that uses a Roe-type Riemann solver, in the context of compressible flow problems) at appropriate polynomial orders, this time with a constraint to enforce a sufficiently strong diffusion for the reflected eigenmodes. DG's numerical diffusion was taken as a reference because DG is regarded as having an adequate balance between dissipation and accuracy. The resulting dispersion and diffusion plots are shown in figure 5 for various orders, with the curves for the spurious modes at the bottom part of the graphs. Once again, the newly designed kernel for the CG-SVV operator, that we named "DG" kernel given its connections to the DG method in terms of diffusion and dispersion characteristics, maintains the Péclet number fixed and guarantees that raising P increases the resolution power per degree of freedom. In order to validate the novel CG-SVV operator based on the "DG" kernel, we used a simple two-dimensional test case that was designed to mimic spatially developing grid turbulence [49,50]. It consists of a rectangular domain with prescribed unsteady flow at the inlet and symmetry (free-slip) conditions on the sides. This domain was meshed with quadrilaterals in a structured grid pattern featuring two regions. The upstream region had square-shaped elements, whereas the downstream region had a four times larger streamwise mesh spacing. The interface between these two regions marked a potentially unstable threshold zone across which mesh spacing changed abruptly. Results obtained without SVV for P = 8 and n = 11 elements in the crossflow direction are shown in figure 6. In particular, the vorticity contour plot obtained without SVV, upon close inspection of the upstream region to the mesh coarsening line, shows the presence of spurious reflections originating at the threshold zone. This case diverged shortly afterwards, along with other test cases where we used different polynomial orders, ranging from P = 3, to P = 7 (while maintaining the same number of degrees of freedom -i.e. decreasing the number of elements). In particular, at lower polynomial orders, higher Reynolds numbers were required for the simulation to diverge, while for higher polynomial orders, the test cases started to diverge at correspondingly lower Reynolds numbers. When the novel DG-based CG-SVV operator was added, however, up to polynomial order P = 8, all cases remained stable no matter how high the Reynolds number tested. Furthermore, no spurious reflections have been noticed. And although solutions at lower orders (P = 3 in particular) might have seemed excessively regularised by the "DG" kernel SVV, superior resolution power was effectively recovered at higher orders, such as P = 4 and above [53]. The current best practice in N ektar has been to use the "DG" kernel in general while employing the "power-law" kernel only along homogeneous directions in case periodic boundary conditions are present. In either case, SVV's base coefficient takes the final form where Pe * is a reference constant Peclét number (typically of order one), v and h are respectively local measures of advection velocity and mesh size (along the relevant directions of flow propagation), and P is the chosen polynomial order. The entries defining the "power-law" kernel are those given in (3.13), whereas those of the "DG" kernel have no explicit formula (as the output of an optimisation procedure), having their values given in table form for each order P in [50]. Furthermore, we note that SVV's kernel operation from (3.11) is accommodated in polynomial space as in whichĝ denotes the vector of element-wise coefficients in the polynomial expansion of ∇u, while R is a "rotation" matrix that projectsĝ onto an orthogonal (Legendre) hierarchical basis, and Q is the diagonal matrix whose entries are given by the chosen SVV kernel function. Therefore, the proposed CG-SVV operator has the final form This CG-SVV operator has been implemented in Nektar++ and extensively tested in different test cases of increasing complexity, always displaying strong robustness and preventing numerical instabilities. This approach offers a decent balance between robustness and accuracy for high-order discretizations, being effectively a parameter-free strategy from the practitioners's perspective and allowing for under-resolved turbulence simulations over complex geometries at high Reynolds numbers. This approach, in conjunction with the dealiasing techniques described in section 3.1, was used to obtain the numerical simulations depicted in section 4. Without the concurrent use of dealiasing and SVV regularization, the numerical simulations presented in this paper would have diverged due to numerical instabilities.

Computational efficiency.
One of the advantages of higher order polynomial expansions is their ability to exploit high arithmetic intensities; that is, a greater number of floating point operations is performed for each byte of data retrieved from memory at higher orders than at lower orders. Although this naturally makes operator evaluations at higher order more expensive per degree of freedom, on modern computing hardware this effect frequently tends to be masked, as the cost of transferring data from memory to the compute unit tends to be far greater than that of performing floating point operations. In essence, this can make floating point operations nearly 'free' compared to the transfer of data.
However, the efficient implementation of the building-block finite element operators that make up the discretised version of the incompressible Navier-Stokes equations (3.1) that were used in this work, such as inner products or derivatives, still requires careful consideration. Even small changes in polynomial order can yield a large effect on performance due to the many dependent factors, including both the choice of hardware and the problem under consideration. Our recent work in this area [58] has focused on the development of techniques that can automatically determine the most efficient kernels for these operators at runtime. These kernels are aimed to target both clear performance factors (such as operator type and polynomial order), as well as more subtle aspects, such as whether the choice of basis and whether mathematical techniques such as sum-factorisation can be leveraged to increase performance, and the grouping of elements with equivalent numerical properties to improve performance evaluation. The use of these kernels therefore allows an efficient on-node implementation across a broad range of polynomial orders.
For parallelisation between nodes, Nektar++ uses a parallelisation model based on distributed MPI. The METIS library [39] is first used to partition the domain based on the dual graph representation of the mesh. Since, for industrial geometries, the mesh generation process yields a hybrid mesh of prismatic and tetrahedral elements, these are weighted appropriately in the partitioning process, so as to approximately equally distribute the work between processors. Both the pressure Poisson equation and the three scalar velocity correction Helmholtz equations that comprise the operator splitting scheme are solved implicitly, leading to the need for efficient solvers for large, distributed linear algebraic positive definite systems. To achieve this, we apply an iterative preconditioned conjugate gradient (PCG) method for each scalar system. We note, however, that the usual approach taken at lower orders of representing this system as a large, distributed sparse matrix is not generally efficient at higher orders. Instead, when the action of the matrix multiplication is required during PCG iterations, we implement this as many smaller, locally dense elemental matrices, combined with a sparse gather operation to apply C 0 connectivity. The gather operation uses the efficient Gslib library in parallel, based on the code Nek5000 [18]. Our tests have demonstrated that this combination allows both performance and scalability of the PCG solver up to 10 5 cores at high levels of efficiency. Finally, we discuss preconditioning of these methods. All four PCG solvers use an effective and scalable low-energy basis preconditioner [73,28] and, optionally, projected initial conditions based on storage from solutions at previous timesteps [20]. The solution of the pressure Poisson equation, however, is particularly ill-conditioned, and therefore requires the use of coarse-space preconditioning, which requires the inversion of the linear finite element representation of the mesh. Currently, the coarse space preconditioner uses a direct parallel inverse method developed by Fischer and Tufo [79].
The combination of all of these intra-and inter-node performance optimisations, together with a robust meshing pipeline and stabilisation strategies, has yielded an efficient solver at high orders that is capable of tackling complex configurations required for industrially relevant problems at superior computational performance.
4. High-order iLES for the aerodynamic (re)design of the Elemental Rp1 road car. The developments in terms of high-order mesh generation and numerics implemented in NekMesh and Nektar++ that were described in sections 2 and 3 were used to perform large-scale high-fidelity LES simulations of a real road car, namely the Rp1 car [16]. The Rp1 car, produced by the British manufacturer Elemental Cars since June 2016, is a high-performance, light-weight, road-legal vehicle that weighs 580 kg. It is equipped with a 320 bhp engine and accelerates from 0 to 60 mph in just 2.6 seconds. Its current configuration, depicted in Figure 7 along with its digital twin, delivers 400kg of downforce at 150 mph. The main targets of the high-fidelity LES simulations performed were to analyze and improve the external aerodynamic performance of the car. In particular, we analyzed three configurations, namely designs 1, 2, and 3, that we denote as D1, D2, and D3, respectively. D1 is the baseline design, while D2 and D3 are two design upgrades that were implemented to improve aerodynamic performance and pilot comfort. This study was a collaborative effort put forward by Imperial College London, Elemental cars group [15], London Computational Solutions (LCS) [42] and Silicon Graphics International Corp (now Hewlett Packard Enterprise). All the simulations were carried out using the open-source library Nektar++. In the following, we describe the configuration of the simulations and we discuss qualitatively the results obtained. 4.1.1. Mesh generation. The hugely complicated nature of the geometry under consideration, which involves hundreds of CAD surfaces and thousands of edges, presented a significant test for the mesh generation process outlined in section 2. The high-order meshes used in this study ranged from approximately 2.2 million to 4 million elements, and they are comprised of a combination of triangular prisms in the near-wall boundary layer region with 6 layers used to resolve separation accurately, and tetrahedra in the rest of the domain. We adopted three different meshing strategies for the three car designs, D1, D2, and D3. In particular, we initially used the first and second analytic curving method described in section 2.3, for the surface mesh generation stage. The first analytic curving method was applied to the first design, namely D1, where the creation of a viable mesh took approximately 4 months of iterations between CAD healing and high-order meshing. The final mesh was composed of 2.2 million elements, including 600000 prisms within the boundary-layer mesh. The second analytic curving method was instead applied to the second design, namely D2, that was constructed based on the results from D1. The mesh for D2 contained 2.6 million elements, including 700000 boundary layer prisms. While, both analytic curving strategies were effective in generating viable meshes for D1 and D2, they were deemed excessively user-driven and time-consuming. The mesh for D3 was obtained with what we named the projection curving method, described in section 2.3. Indeed, this method, along with the use of CADfix for geometry repair, Star-CCM+ for linear and boundary-layer mesh generation and NekMesh for high-order mesh generation was the least time-consuming and user-driven. On the one hand, it allowed cutting the overall meshing time from 4 months to 1 week. On the other hand, the mesh was created almost without user effort, and the process required one single execution of linear mesh generation via Star-CCM+ and one single execution of NekMesh. No repetitive cycles of CAD healing were required. The final mesh produced for D3 contained approximately 3.5 million elements, including 1.24 million boundary-layer prisms. Although the meshing strategy for D3 is the most robust and efficient, we should point out that allows for inaccuracies in the geometry description, as a small percentage of surface elements are included in the "exclusion set", thereby remaining straight-sided. However, in the results that we present in section 4.2.2, no obvious unphysical results could be noticed, and the flow physics in those regions of the car that were geometrically identical between D2 and D3 was indistinguishable. Figure 8 shows the three surface meshes for D1, D2, and D3 obtained via the three different strategies just outlined. It is possible to note that the surface mesh for D2 and D3 has a smoother element distribution than D1, due to the enhanced triangulation strategies adopted, namely the second analytic curving, and projection curving. We remark that for all the simulations carried out we did not include the spokes that are instead shown in figure 4.
(a) (b) Fig. 8. Surface meshes of the three car designs considered, D1, D2, and D3. Figure (a) shows a comparison between D1 and D2, while figure (b) shows a comparison between D1 and D3. It is possible to notice a smoother distribution of surface elements for both D2 and D3 when compared to D1. This is due to the different approaches adopted for the surface mesh generation, where for D1 we used the first analytic curving, while for D2 and D3 we used the second analytic curving and projection curving, respectively.

Simulations setup.
The simulation setup required to analyze the aerodynamic of the Rp1 car was challenging due to the high Reynolds number, Re = 10 6 , based on the length of the car L, and to the complex geometries involved. In order to obtain stable yet accurate simulations, we adopted some simplifications that are commonly adopted when performing the aerodynamic car simulations.
The Rp1's unique double diffuser design and the the road-tyre interaction were carefully modeled with some simplifications that allowed a significantly easier simulation setting while not undermining the accuracy of the results. Rp1's double diffuser design has intakes and outlets on the body of the car. To model these accurately, we extracted the corresponding surfaces and constructed an appropriate velocity boundary field, by enforcing a prescribed average normal velocity on each surface. The boundary condition has to also satisfy the no-slip condition at the edges of the inlet surfaces and this was imposed by solving a Helmholtz problem where the positive Helmholtz coefficient was chosen to provide a layer of desired thickness at the edges. This was possible by using a manifold solver developed by Cantwell et al. [6]. In addition, the road-tyre interaction was modeled using a contact patch along with rotational boundary conditions to a simplified tyre geometry. This allowed efficiently modeling the rotation of the tyre without the need for a moving mesh. Finally, we simulated only half of the car, by using a symmetry plane. These simplifications allowed a faster simulation setup and simulation runtime, without altering the flow features that were of interest in this study.
The simulations were initialized adopting a similar strategy to Lombard et al. [43], where we began at a lower Reynolds number, Re = 10 4 . The initial conditions were imposed impulsively so that u/U = 1, where U is a reference velocity, and all nonmoving surfaces are treated as no-slip boundaries. At the outflow, we imposed the boundary condition developed by Dong et al. [14], which balances the kinetic energy influx through the outflow boundary condition to prevent instability. The Reynolds number is then gradually increased by a factor of ten until the target Reynolds number, Re = 10 6 , is reached. Between increases, the flow is evolved for two convective time units t c = L/U to allow the flow to adapt to the new flow physics dictated by the higher Reynolds number. We note that the Reynolds number used in this study, Re = 10 6 , is still lower than the experimental value of approximately Re = 2 × 10 6 . Nevertheless, to the authors' knowledge, this is the highest Reynolds number ever attempted for a real car geometry, using iLES in conjunction with spectral/hp element methods. The simulations were run at a polynomial order 4 (fifth order of accuracy) with a second order implicit-explicit time-stepping. To increase the accuracy of the integration, a seventh order of quadrature was used, that consisted of the so-called global consistent integration introduced in Mengaldo et al. [48], and described in section 3.1. This allowed us to mitigate aliasing-driven instabilities. We note that the use of a higher quadrature meant that the meshes were generated at a higher order of P = 7 to ensure the elements were valid with this quadrature. In addition, the SVV kernel presented in section 3.2 was used to filter the under-resolved scales and enhance the stability of the simulations. Without the use of both, consistent integration and SVV was not possible to obtain numerically stable calculations.
In terms of computational requirements, the simulation were performed in parallel using 1000 to 4000 cores, depending on node availability on a SGI UV machine, and it took approximately 3 days to evolve the flow for one convective length. The high-fidelity flow simulations performed allowed unique insights into the aerodynamic performance of the Rp1's complex geometry, highlighting a number of important regions for aerodynamic optimisation. These results are discussed next.

4.2.
Results. In this section, we discuss the results obtained on the three car designs introduced earlier, D1, D2 and D3. In particular, in section 4.2.1, we present a comparison between D1 and D2, while in section 4.2.2, we present a comparison between D1 and D3. The results focus on iso-contour quantities that allowed us to identify important flow features suitable to optimization. 4.2.1. First design stage. The first simulation performed was for design 1 (D1), that is the first baseline model of the car that was produced, followed by the simulation of design 2 (D2). The latter design, D2, implements an upgrade of the aerodynamic package with respect to D1. Figure 9 shows the meshes and associated results for D1 and D2. In particular, figure 9-(a), shows the meshes for the two designs (that are the same as depicted in figure 8), while figure 9-(b) shows the associated flow results. Figure 9-(c), shows the same flow results as the right image in figure 9-(a), from a sideview perspective. The flow results depicted in figure 9 show iso-contours of the total pressure coefficient Cp 0 = 1 − (u/u ∞ ) 2 = 0, and they are colored as pressure. From the high-order D1 simulation two key findings were made, neither of which were identified in low-order RANS simulations [77]. Firstly, there appeared to be strong vortical structures hitting the drivers helmet (see figures 9-(b) and 9-(c) for D1). Indeed, aerodynamic disturbances were also reported by the test drivers, confirming the finding of our iLES. These vortical structures were generated by the windscreen, and the issue was fixed in D2 with a redesigned console area in front of (a) (b) (c) the driver, significantly enhanced the driving experience. The D2 iso-contour shows clearly that these structures now pass cleanly over the drivers head and are less noisy in terms of small-scale vortices (figures 9-(b) and 9-(c), for D2). Secondly, the roll hoop produced significantly more drag and separation than predicted. D2 had a redesigned roll hoop with an aerofoil profile as opposed to a circular cylinder. This aerofoil profile was slightly angled to help control the flow over the new Gurney flap at the rear of the car. This, combined with redesigned diffusers on the underside, led to increased downforce over D1. The trend of changes in downforce and drag between the D1 and D2 simulations was well predicted by the high-order simulations and matched the trend in the RANS results [77].

Second design stage.
The final simulation performed in this study is for design 3 (D3), which is a full aerodynamic upgrade over the previous two designs, D1 and D2. This car is designed to achieve extremely high levels of downforce specifically for track racing. This design includes a fully redesigned floor, front splitter and the addition of a rear wing. The ride height has also been altered, raising the car at the rear for increased diffuser performance and lowering the front of the car to increase in-ground effect of the front splitter. Figure 10, similarly to figure 9, shows again the total pressure coefficient, colored by pressure. The results are consistent with the D2 simulation and are shown alongside the D1 results. The trends in increasing downforce were well predicted by the highorder simulation. In figure 10, we note that it is possible to observe a noticeable offset in the geometries of the two cars, due to the alteration in ride height.
Although we used an approach for meshing D3 that can be aggressive in preserving mesh validity versus geometrical accuracy, thereby leaving a non-negligible number of elements straight-sided as opposed to curvilinear (the number of straight-side elements was < 5%), the flow physics is in agreement with D1 and D2, and no unphysical results could be readily noticed. However a more comprehensive study should be undertaken in order to draw stronger conclusions on whether the compromised geometric accuracy of our most robust meshing pipeline is a compromise worth taking. The early results here show that it may well be.

5.
Conclusions. This work demonstrates the significant progress in leveraging Nektar++ for the goal of moving high-order CFD from academia to industry, and some of the benefits that high-order methods can provide to the industrial community. Significant developments had to be made in order to overcome obstacles in both mesh generation and solver technologies.
The novel high-order meshing workflow described in section 2 allowed the robust, efficient and accurate high-order mesh generation of an extremely complex geometry. This was a key to prevent unphysical diffusion being generated at no-slip surfaces that constitute the geometry of the car, and to avoid using an excessive number of linear elements in proximity of curved geometrical features. Indeed, the meshing framework has proven reliable and viable for industrial geometries, and it was a key enabler to successfully accomplish the simulations presented in section 4.
Along with the high-order meshing framework, the novel numerical technologies described in this section 3 were essential to successfully produce the CFD workflow described in this paper. More specifically, the consistent integration for mitigating aliasing-driven instabilities in conjunction with the CG-SVV based on the "DG" kernel for treating under-resolved scales were the keys to obtain stable simulations and accurate simulations. The use of the novel CG-SVV approach, in particular, was instrumental to produce high-fidelity results. In fact, the CG-SVV operator is based on ESA's estimates for spectral/hp element methods, where numerical dispersion and diffusion characteristics are quantified in wavenumber/frequency space, as in classical von Neumann analyses. Although mostly feasible for simple model problems, such as one-dimensional linear advection, ESA's estimates for spectral/hp methods have been shown to hold more generally for nonlinear problems and even turbulence simulations [55,53,17]. These estimates quantify how the numerical error (diffusion in particular) increases at large wavenumbers/frequencies, affecting especially the smallest, poorly-and under-resolved scales. When interpreted alongside implicit LES experiments, ESA reveals that spectral/hp discretisations of higher polynomial orders perform best because they affect the large and intermediate flow scales as little as possible, while gradually introducing a sufficiently strong diffusion to regularise the smaller scales. This way, the governing equations and thus the turbulence physics are solved more accurately for a wider range of scales, rather than having to rely on (often restrictive) modelling hypothesis. We note that past studies also focussed on how numerical diffusion acts as a turbulence model [45,56,9,10,52]. However the ESA framework was deemed the most promising route to achieve high-fidelity implicit LES results.
The simulations in section 4 represent the first high-order (having solution polynomials = 4) accurate results of a real automotive geometry by means of implicit LES technologies using spectral/hp element methods. They represent one of the first attempts at bridging the gap between applied mathematics academic research and industrial applications within the context of CFD.