Industry-Relevant Implicit Large-Eddy Simulation of a High-Performance Road Car via Spectral/hp Element Methods
Gianmarco Mengaldo, David Moxey, Michael Turner, Rodrigo C. Moura, Ayad Jassim, Mark Taylor, Joaquim Peiro, Spencer J. Sherwin
IINDUSTRY-RELEVANT IMPLICIT LARGE-EDDY SIMULATION OFA HIGH-PERFORMANCE ROAD CAR VIA SPECTRAL/ HP ELEMENT METHODS ∗ GIANMARCO MENGALDO † , DAVID MOXEY ‡ , MICHAEL TURNER § , RODRIGO C.MOURA ¶ , AYAD JASSIM (cid:107) , MARK TAYLOR , JOAQUIM PEIR ´O § , AND
SPENCER J.SHERWIN § Abstract.
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 charac-terized by high Reynolds numbers and complex geometries. In particular, we describe the numericalmethods, software development and steps that were required to perform the implicit LES of a realautomotive car, namely the Elemental Rp1 model. To the best of the authors’ knowledge, this simu-lation 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 envelopecurrently allowed in industry, where steady-state modelling remains the standard. To this end, anumber of novel developments had to be made in order to overcome obstacles in mesh generationand solver technology to achieve this simulation, which we detail in this paper. The main objectiveis to present to the industrial and applied mathematics community, a viable pathway to translateacademic developments into industrial tools, that can substantially advance the analysis and designcapabilities of high-end engineering stakeholders. The novel developments and results were achievedusing the academic-driven open-source framework Nektar++ . Key words. spectral/ hp element methods; high-order CFD; high-order mesh generation; CADintegration; high Reynolds number flows; high-fidelity simulations; implicit LES; under-resolvedDNS; industrial applications; automotive design. AMS subject classifications.
1. Introduction.
The accurate modelling and prediction of fluid dynamics struc-tures and aerodynamic forces around complex bodies is a classical and challengingproblem that lies at the heart of modern automotive design. The increasing availabil-ity and lowering cost of computational modelling compared to experimental testingmeans that computational fluid dynamics (CFD) now forms a key part of the engi-neering design process. When considering the high Reynolds number flows aroundthe complex bodies found in these problems, there are two key considerations thatneed to be taken into account: (a) creating highly accurate meshes for the geometryat hand, and (b) producing physically-faithful results, while using numerical discreti-sations that are necessarily unable to capture all the fluid dynamics scales involveddue to the large computational costs required.In industry, the most common practice revolved around the use of ensemble-averaged steady-state modelling, commonly using either Reynolds-averaged Navier-Stokes (RANS) or, more recently, Detached eddy Simulation (DES). Typically these ∗ Submitted to the editors DATE.
Funding:
EPSRC under the Platform Grant PRISM: Underpinning Technologies for FiniteElement Simulation (EP/L000407/1); EU Horizon 2020 project ExaFLOW (grant 671571). EPSRCand Airbus under an Industrial CASE studentship. † National University of Singapore (NUS), Singapore ([email protected]) ‡ University of Exeter, UK ([email protected]) § Imperial College London, UK ([email protected], [email protected],[email protected]) ¶ Instituto Tecnol´ogico de Aeron´autica, Brazil ([email protected]) (cid:107)
Hewlett Packard Enterprise, UK ([email protected])
London Computational Solutions, UK ([email protected])1 a r X i v : . [ c s . C E ] S e p G. MENGALDO, ET AL. are used in conjunction with low-order numerical discretisations, predominantly thesecond-order accurate finite volume method, together with linear meshes of planar-faced elements. Although this approach is computationally cheap, it can underminethe accuracy of the simulations and keep the practitioner from identifying aerody-namic design solutions that would be otherwise feasible if the CFD methodology wasable to more faithfully describe the flow physics [74]. Indeed, on one hand, the useof highly diffusive numerical schemes and steady-state modeling strategies leads tosuppress under-resolved scales and their feedback on the mean flow, thereby addi-tionally corrupting the accuracy of the simulation. On the other, the use of linearmeshes requires high degrees of refinement in order to capture highly-curved surfacesin the geometry, which may introduce artificial diffusion and significantly alter theflow evolution that one would otherwise observe. Addressing point (a) and (b) withinthe context of industrial applications is therefore of primary importance to advancethe state-of-the-art of CFD, and can dramatically enhance the computational aero-dynamic analysis and design capabilities of high-end engineering stakeholders.Point (a) requires substantial theoretical and algorithmic advances in high-ordermesh generation. From this perspective, the aim is to produce a methodology forgenerating high-order meshes that achieve levels of robustness over a range of verycomplex cases comparable to that found in commercial linear mesh generators. In-deed, commercial packages, such as Pointwise [65] and Centaur [8], have attempted toadd 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 in-troducing the geometric curvature into elements lying on the surface. Several researchgroups have been working on such an approach, see for example [13, 73, 44, 71, 69, 84],which gives an indication towards the level of interest in this area. The challenge ofthis approach, however, lies in the second step of the process: introducing curvaturepurely at the surface will typically lead to elements becoming inverted and unsuitablefor computation. This has proven a significant stumbling block in the use of high-ordermethods in industry, and many of the aforementioned works focus on developing tch-niques to regularize or heal the mesh as part of the a posteriori process. While this isthe main route to high-order mesh generation, some alternatives exist. These includethe use of the Nash theorem to map a Riemannian space onto an isometric Euclid-ean space of higher dimension [61, 86, 11], semi-structured approaches to block meshgeneration [3], anisotropic high-order boundary-layer mesh generation [60], and PDE-based approaches based on elastic analogy [67]. In this work, we tackle the problem ofhigh-order mesh generation by using a variational framework that uses the minimiza-tion of an energy functional. The process begins by defining a linear surface mesh thatconforms with the underlying geometry, typically defined in a computer-aided design(CAD) format that follows standard ISO 10303 (informally known as standard for theexchange of product model data, or STEP) [68]. The linear surface mesh is then com-plemented 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 theunderlying geometry, and the interior mesh entities (edges and surface) of the linearvolume mesh are deformed to accommodate the curved boundaries imposed by thegeometry. The latter step is achieved by the minimization of the energy functionalmentioned above, that can be arbitrarily selected by the practitioner depending on themesh requirements. The linear mesh is obtained beforehand, and uses a combinationof the commercial packages
CADfix and
Star-CCM+ for CAD healing and for linearand boundary-layer mesh generation, respectively, as well as
NekMesh [79], which is
NDUSTRY-RELEVANT IMPLICIT LARGE-EDDY SIMULATION
Nektar++ project [62, 5, 58].The second point, (b), is instead closely linked to the impossibility of carryingout 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 powerrequired to perform a fully resolved flow simulation, also known as direct numericalsimulation (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 acomputer can perform as an approximate measure of computational power, and wecan then model the time to solution under different regimes of computational powergrowths. This gives a rough estimate of the number of years required to achieve DNScapabilities for different Reynolds numbers. Figure 1 depicts two different computa-tional power growths. The red line assumes doubling flops every three years, while theblue line assumes a doubling of FLOPS every one and a half years. They represent aless and more optimistic scenario than Moore’s law, respectively. We can appreciatehow, under the most optimistic power growth scenario (blue line), DNS capabilitiesfor automotive simulations, which corresponds to Reynolds numbers of the order of10 , will be reached approximately by 2040. Figure 1 is an ideal scenario, as it doesnot take into account additional costs incurred, and it also assumes using the topsupercomputers, as soon as they become available. The latter aspect is obviously notachievable for analysis and design purposes in an engineering workflow. Therefore, theability to perform DNS for industrial purposes at high Reynolds numbers on a regularbasis is, at least, two or more decades ahead, even with current advances in compu-tational technologies (e.g. many-core architectures, graphical processing units, andquantum computing). Given these limitations, we need to deal with under-resolutionfor the next several years, while pushing the boundaries of high-Reynolds number flowphysics for industrial problems. Fig. 1 . Computational resources required to perform the direct numerical simulation (DNS) asa function of the Reynolds number [37]. Yellow and red lines represent projections on when DNScapabilities are due to become available as computational power increases following two differentMoore law’s paths.
A significant body of work has been produced in the past several years to deal
G. MENGALDO, ET AL. with high-Reynolds number problems. As briefly mentioned above, RANS and DESare frequently used in the context of industrial CFD, see for example [23, 85, 4], butthey typically fail to capture unsteady flow features, which may be crucial for theanalysis and design of a given aerodynamic solution. This aspect is especially markedwhen these two modeling choices are adopted in conjunction with low-order numericalschemes, that are commonly characterized by significant diffusion and dispersion,thereby causing an unphysical behaviour of the flow physics [7]. In order to overcomethese issues, several research groups have been working on LES, both physical (orexplicit) and implicit, over the past few decades. Explicit LES uses physical principlessuch as eddy viscosity to devise closure models that can approximate the under-resolved scales, while implicit LES utilizes the diffusion introduced by the numericaldiscretization. For a detailed review on LES methods, the interested reader can referto Pope [66]. LES is particularly promising for automotive applications, since theseregimes heavily feature flow physics that steady-state approaches struggle to capture,such as vortex separation and interaction [74]. In addition, when used in conjunctionwith low-diffusive and dispersive numerical scheme, the resolution power is enhancedand can produce unprecedentedly accurate flow simulations. In this work, we use animplicit LES approach, where we treat the under-resolved scales via a novel artificialnumerical regularization approach based on a so-called spectral vanishing viscosity(SVV) kernel, and we adopt robust dealiasing techniques to mitigate aliasing drivennumerical instabilities that are normally present due to the low-diffusive behavior ofthe numerical framework. In particular, we show that the use of SVV and dealiasingstrategies, along with key improvements in the computational efficiency of kernelcalculations, 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 ofthe Rp1 automotive car at a Reynolds number of 10 based on the car length. Allthe simulations performed use the incompressible Navier-Stokes solver implementedin the open-source spectral/ hp element framework Nektar++ .In this paper, we present the key advancement in high-order meshing and simula-tion technologies that we implemented in
NekMesh and
Nektar++ in order to achievethe challenging simulation of the Rp1 car, with the aim translating academic researchefforts to industrial CFD applications. In particular, section 2 describes the high-ordermeshing strategy we adopted, focusing on the steps required for surface and linearmesh generation as well as on the a posteriori high-order deformation. In section 3,we describe the novel simulation technologies implemented in terms of regularisationapproaches for the under-resolved scales, dealiasing techniques for numerical stability,and computational efficiency. Section 4 presents the high-fidelity simulation of theRp1 car, and finally, in section 5, we draw some final remarks.
2. High-order mesh generation for complex geometries.
The startingpoint 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 andhigh-quality mesh. Mesh generation is a complex and involved process, particularlywhen the geometry under consideration is a complex, three-dimensional body, sincethe important regions within the domain to resolve from a fluid dynamics perspectiveare often unknown a priori . The process may therefore require many iterations ofmeshing, simulation and analysis in order to arrive at a final mesh deemed suitableenough for accurate results. Mesh generation tools are therefore required to be bothflexible and robust, leading to a wide variety of approaches in industry to deal with
NDUSTRY-RELEVANT IMPLICIT LARGE-EDDY SIMULATION a posteriori mesh generation.
Our frame-work to produce a viable high-order mesh, similarly to other state-of-the-art methods[13, 73, 69, 32, 84], is based on an a posteriori approach. In this setting, one firstgenerates a coarse, linear mesh of the desired geometry. From this, high-order nodesare added within each element. Any nodes that lie on curved surfaces are projectedonto 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 typ-ically adopt an isoparametric approach, in which the same basis functions that areused to represent high-order solutions within each element are also used to define thegeometric coordinates within the element. Such a mapping can be expressed thereforeas the summation(2.1) x = φ M ( ξ ) = N (cid:88) n =1 x n (cid:96) n ( ξ ) , where x is the spatial location of the a point within an element, ξ is a corresondinglocation within a reference element, x n is the location of a given node in the element,and (cid:96) n is the polynomial interpolant (e.g. a Lagrange polynomial interpolant). Inparticular, mesh entities, namely faces and edges, on the boundary of the geometryare curved by adding high-order nodes, and then enforcing that the nodes lie on theboundary. Moreover, the position of the nodes x n must be adjusted so that theresulting mapping is a faithful polynomial representation of the geometry [73].Although from this description this may seem a simple task, for complex geome-tries that can have several hundred curved surfaces, the ability to create valid high-order elements is extremely challenging. In many cases, when simply curving the edgeand/or face in proximity to a given curved geometrical feature, it is relatively easy togenerate invalid elements. Figure 2 shows seven different mesh elements in proximityof a curved surface, depicted in grey. The first mesh element (a) is a valid linear ele-ment that completely misses the curvature of the underlying surface. Indeed, in orderto capture the curvature of the surface with a linear mesh it is necessary to add more G. MENGALDO, ET AL. elements around the surface than would be necessary with high-order elements; thelatter benefits from a superior accuracy that allows for larger element sizes than linearelements. The second element (b) is a valid high-order element, where only the edgein proximity to the geometry is curved. The third and fourth elements, (c) and (d),are both invalid high-order elements, due to near tangent (c) or intersecting edges (d).Mesh elements (b), (c), and (d), clearly highlight the challenge of high-order meshgeneration. The curvature of the geometry dictates the size of the element one canhave, and in several occasions, it could lead to invalid elements. Mesh elements (e)and (f) show two possible solutions to the invalid elements (c) and (d). Mesh element(e) deforms the interior edges to accommodate the curvature of the boundary, while(f) uses a quadrilateral element to increase the permitted curvature. In both cases,there are no longer near tangent or intersecting edges. From the examples in figure 2,
Curved geometry
Valid high-order element
Invalid high-order element (near tangent edges)
Invalid high-order element (intersecting)
Valid high-order element via interior deformation
Valid high-order element via increased permitted curvature (a)
Valid linear element (b) (c) (d) (e) (e)
Fig. 2 . Different valid and invalid high-order elements. it should be immediately clear that the challenge of high-order mesh generation ismultifold. On the one hand, it is necessary to create curved mesh entities that con-form to the underlying curved geometrical features (where the geometry is typicallyprovided in a CAD format). On the other hand, it is essential to develop a strategyto handle invalid high-order elements and produce a valid high-quality mesh.Our high-order meshing strategy handles the challenges just described and con-sists of three steps, (1) the definition of the geometry in such a way that is canbe efficiently queried for curvilinear mesh generation, (2) the generation of a high-quality linear mesh, and (3) the generation of curved elements. These three crucialsteps are described in detail in sections 2.2, 2.3 and 2.4, respectively. The correspond-ing meshing framework has been implemented in the code
NekMesh , that is part ofthe open-source
Nektar++ library [5]. For a more comprehensive explanation of themesh generation process, the interested reader can refer to Turner et al. [80].
The first stage of our meshing approach is the def-inition of the geometry in a way that ensures its accurate representation and permitsgeometrical enquiries during the various mesh-generation stages. The most commonapproach to define an arbitrary complex geometry is by a boundary representation(also referred to as “BRep”). The boundary representation contains surface patchesand curves that link these surface patches. The surface patches and curves defin-ing the underlying geometry are typically topologically connected, thereby forminga closed shell. These surfaces and curves are parametrised onto a lower-dimensionalspace, that is the one-dimensional space x ( s ) for curves, and the two-dimensionalspace x ( s , s ) for surfaces, where x = ( x, y, z ) denotes the three-dimensional coordi-nate of a point in the boundary representation of the geometry. The main geometricalenquiries required for mesh generation consists of three key operations: NDUSTRY-RELEVANT IMPLICIT LARGE-EDDY SIMULATION • parametric coordinate to real space location, i.e. determining x from ( s , s ); • projection of a location in space onto a boundary representation entity, i.e.determining ( s , s ) from x ; and • calculation of the first- and second-order derivatives of the parametric map-ping with respect to s and s .In addition, in this work we do not consider the effects of incomplete or low-qualitygeometric representations; The boundary representation must be watertight and shouldrepresent accurately the geometry under investigation. These two requirements arecrucial for high-order mesh generation, as they impact the number of high-order validelements 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 stan-dard CAD (computer-aided design) format that supports a boundary representationdescription, such as STEP. In several CFD applications, this assumption is fairly com-mon, especially among industry practitioners. The CAD object is then probed via NekMesh , that performs geometrical queries through interfaces to a CAD engine. Togive support for different CAD engines and implementations,
NekMesh implements alightweight interface to the key operations above. Concrete implementations of thisuse CADfix [34] and OpenCASCADE [63], both of which permit reading CAD datain various formats, including STEP. In particular, the checks and queries on the CADobject were performed by
NekMesh , while the correction of any errors was performedby CADfix which provides an extensive toolset for fixing geometrical inconsistenciesin CAD data. CADfix is a commercial software, that provides extensive CAD han-dling 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 usefulfor decomposing the domain into partitions. This decomposition has proven usefulfor designing high-quality boundary-layer meshes near the wall surfaces, aspect thatwas critical for the high-fidelity simulations in section 4. OpenCASCADE is an open-source code with reduced capabilities than CADfix, and represents the default CADhandler within
NekMesh , due to its free-availability, that is at the core foundations ofthe
NekMesh and
Nektar++ projects. Nevertheless, the use of CADfix in this workwas essential to obtain a robust and efficient high-order meshing pipeline, and ulti-mately obtaining valid high-quality meshes. Access to CADfix, or some alternativeapproach to geometry healing and defeaturing, is a fundamental requirement for thedevelopment of a robust and successful meshing strategy.Once the CAD representation of the geometry is satisfactory, it is possible tocomplete the definition of the computational domain, including the external bound-aries, and proceed to the next meshing stage, that is generating a valid linear mesh.This step is described next.
The generation of a high-quality linear mesh from the CADrepresentation of the geometry is a demanding task. Indeed, it is necessary to have asufficiently coarse initial mesh that serves as a starting point for the high-order mesh-ing process. This initial mesh is typically obtained via a user-driven manipulation ofthe mesh spacing definition at all points within the computational domain to achievesatisfying geometric accuracy, coarseness and mesh gradation. This procedure is typi-cally time-consuming and challenging, possibly involving a large number of iterationsbefore 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 advancedlinear meshing capabilities of the commercial package
Star-CCM+ with the high-ordermeshing procedures of
NekMesh . This allowed us to produce an automatic industrial
G. MENGALDO, ET AL. linear meshing pipeline, that required minimal intervention from the user. The mainobstacle that we encountered when using
Star-CCM+ , and that is common to encounterwhen using commercial CFD software, was the lack of information relating the linearmesh to the original CAD definition. Often the CAD surfaces are replaced by a finetriangulation that represents the “true” boundary of the geometry within a giventolerance. The linear mesh is then built using this approximate representation andthe original CAD information discarded. This is a significant issue since inaccuraciesin the definition of the geometry could lead to errors in the high-order mesh generationprocess, 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 theCAD surfaces to reconstruct the missing CAD information. In particular, we importedthe linear mesh generated by
Star-CCM+ into
NekMesh . The import was performedby incorporating curvature-driven mesh size restrictions such that the boundary ofthe 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 isprocessed in order to obtain parametric coordinates in the CAD representation thatwill be used to enhance the adherence of the linear mesh to the CAD model. Fromthis 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, theCAD model itself, and the “auxiliary triangulation”. The second step is to definea bounding box (inflated by 5% in each direction), for each CAD surface, and storethese boxes into a k -dimensional ( k − d ) tree data structure [1]. Once the two CADrepresentations 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 ofpotential parents surfaces. The use the k − d tree makes this process particularlyefficient and reduces the number of potential parent surfaces to just a few candidates,typically in the order of 2 or 3. The surface mesh node is then projected onto the CADsurface for each of the parent surface candidates. This is achieved by identifying thenearest vertex in the “auxiliary triangulation” of the CAD surface, and then using thisvertex 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 tothe projected node.We note that, since the “auxiliary triangulation” of the CAD surfaces is not di-rectly related to the Star-CCM+ linear mesh, the surface node of the linear mesh mightnot lie on the CAD surface. Therefore, once the parent CAD surface associated tothe 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 therequired displacement of the node is greater than 10% of the length of the edges inthe local mesh, and if the movement of the node induces inverted (i.e. invalid) ele-ments. In these two instances, the node is added to an “exclusion set” with a view topreventing deterioration of mesh quality or generation of invalid elements. The finalproduct of the projection process just described is a set of nodes that are accuratelyplaced 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 allthe parametric information required to generate the final high-order curving of thesurface mesh, with the exception that any mesh entity (edge or face) that has a node
NDUSTRY-RELEVANT IMPLICIT LARGE-EDDY SIMULATION projection curving [79], as it involves the projection of the surface mesh generated via
Star-CCM+ onto an “auxiliary triangulation” of the CAD surfaces. This approach isthe one used in the final industrial meshing pipeline, that was adopted for the lastdesign stage, D3, of the road car in section 4. A boundary-layer mesh was also cre-ated by
Star-CCM+ , and consisted of a prism layer in proximity of the surfaces of thegeometry.We finally point out that we initially attempted two additional alternatives, whichwere 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 andD2, of the road car presented in section 4. In both the analytic curving strategies, allthe parametric information associated with the linear surface mesh was known, hencehigh-order curving of the surface was a relatively simple task, unlike in the projectioncurving approach. The first analytic curving strategy used the linear meshing capabil-ities of
NekMesh to create the surface mesh that was later exported to
Star-CCM+ forthe generation of the volume mesh, and of the boundary layer mesh (constituted by anear-wall macro prism layer). This approach was however extremely time-consumingas it involved dozens of iterations between
NekMesh , Star-CCM+ , and the original CADmodel. The second analytic curving strategy tried to address the bottlenecks of thefirst by allowing the
Star-CCM+ linear mesh generator to create the surface mesh froma linearised CAD surface obtained using
NekMesh . The surface mesh was then forcedto explicitly obey to the boundaries of the CAD surfaces. This is critical in an in-dustrial meshing pipeline, as it is typically impractical to clean complex CAD modelswith thousands of surfaces, so as to guarantee a clean enough CAD to enforce thesurface mesh to obey to the CAD surfaces.
The third stage is the generation ofthe curved surfaces, where the coarse linear mesh generated in the previous step anddescribed 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 meshentities that conform to curved geometric features can be handled via two differentstrategies, projection curving and analytic curving , and we found the former to be theoptimal approach for complex industrial CAD models. While curving the boundariesattached to the geometry is a relatively simple task once the parametric coordinatesof the surface mesh are known and accurate, the creation of valid high-order elementsremains a challenge, as depicted in figure 2. Indeed, one of the key research themesin a posteriori high-order mesh generation is to find suitable deformations of theinterior mesh entities. The ability to deform interior mesh entities allows for bettermesh quality and for a significant reduction of invalid elements, see figure 2-(e). Toobtain suitable deformation of interior mesh entities, there exist a certain number ofalternatives, including methods based on partial differential equations [64, 84, 59, 22,30] and on the minimization of an energy functional [77, 24]. In this work we pursuedthe latter, and developed a variational framework able to encapsulate arbitrary energyfunctionals 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 asdescribed in section 2.3. Following [79], we consider a linear mesh with N el straight-sided elements, that is Ω I = (cid:83) N el e =1 Ω eI , where each element is provided with a high-order polynomial basis function based on standard Lagrange polynomials. The map-ping between the straight-sided mesh, Ω I , and the curvilinear mesh, Ω, is denoted by0 G. MENGALDO, ET AL. ⇠ = ( ⇠ ,⇠ ) ⌦ st standard element ⇠ ⇠ ⇠ n y = ( y , y ) ⌦ eI ideal element y n x = ( x , x ) ⌦ e curvilinear element x n I M Figure 1.1.
A triangular element is used for illustration purposes, but the notationis general and applicable to other element types. On the left the mapof a standard (reference) element ⌦ st onto the straight-sided element ⌦ eI through the mapping I :⌦ st ! ⌦ eI and onto the curvilinear element e :⌦ st ! ⌦ e . The deformation mapping :⌦ eI ! ⌦ e is then definedthrough the composition = eM I . of the domain are curved by adding high-order nodes and then enforcing that thenodes lie on the domain boundary. This simple idea for high-order mesh generationhas a key flaw: the curving process can, and almost certainly will, create elementswhich self intersect or have near tangent vertices in all but the most simple geometricand mesh configurations. Figure 1.2 illustrates this. It shows a valid linear elementin figure (a), a valid high-order element in figure (b), a invalid high-order elementdue to near tangent vertices (c) and a invalid element due to self-intersection (d),based on di↵erent profiles of the geometric edge shown with the red line. The blackedges of the elements can be considered to be interior mesh entities and are thereforenot curved.In this example the problem could easily be solved by curving the other twoedges of the triangle. But whereas the domain definition, usually from CAD, defineshow to curve the boundary mesh entities, no such guidance is available on how tocurve the interior entities. Therefore the task of high-order meshing is to obtaina deformation such that the resulting curvilinear mesh is geometrically conformingwith the boundary and its interior is curved to accommodate the boundary defor-mation and produce a valid mesh. However, it is notoriously di cult to obtain a Fig. 3 . Deformation of a straight-sided triangular element, referred to as standard element,onto an ideal element via mapping φ I : Ω st → Ω eI , and onto a curvilinear element via mapping φ M : Ω st → Ω e . The mapping of the deformation φ from the ideal element and the curvilinearelement is defined as the composition of the two mappings, φ : φ M ◦ φ − I . Figure taken from Turneret al. [80] φ : Ω I → Ω.With reference to figure 3, we denote the coordinates of the standard (also referredto as reference) element as ξ ∈ Ω st , those of the curvilinear element as x ∈ Ω e , andthose of the ideal element as y ∈ Ω eI . The mappings between the elements in figure 3are constructed in an isoparametric fashion, such that the nodes ξ n defining theLagrange 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 bereadily generalized to any element shape, in both two and three dimensions.Once the nomenclature for the mesh is specified, we can define the energy func-tional adopted throughout this work, that is(2.2) E ( ∇ φ ) = (cid:90) Ω I W ( ∇ φ )d y , where W depends on the deformation gradient tensor ∇ φ ( y ) = ∂ φ ∂ y and its determi-nant J = det[ ∇ φ ], also referred to as the Jacobian. The functional W can assume var-ious forms depending on the associated energy adopted. We have implemented func-tionals associated to linear elasticity energy, isotropic hyperelasticity energy, Winslowequation energy, and a distortion-based energy. Details of these functionals can befound in Turner et al. [80] and in Turner [79], and are not reported here for the sakeof brevity.Since the ideal mesh is a union of many elements, the functional can subsequentlybe defined element by element(2.3) E ( ∇ φ ) = N el (cid:88) e =1 (cid:90) Ω eI W ( ∇ φ )d y = N el (cid:88) e =1 W [ ∇ φ M ( ξ ) ∇ φ − I ( φ I ( ξ ))] det( ∇ φ I )d ξ , where we used the composition φ = φ M ◦ φ − I . To fully define equation (2.3), weneed to specify the expression of φ M and φ I . The first, φ M , assumes the usual formof an isoparametric mapping commonly found in nodal spectral element methods andspecified in equation (2.1). The second, φ I , assumes an analytical form that is a NDUSTRY-RELEVANT IMPLICIT LARGE-EDDY SIMULATION v , v , and v , the mappingis given by(2.4) φ I = ( v − v ) ξ + ( v − v ) ξ + v . The minimization of the functional is then formally defined as follows(2.5) find min φ E ( φ ) = min φ (cid:90) Ω I W ( ∇ φ )d y , where we seek the configuration of nodes of the linear mesh, denoted by N =[ n , . . . , n N nodes ], that minimizes the energy functional in equation (2.5). In prac-tice, this means to start from the initial distribution of nodes in the linear mesh, N , in which the curvature of the boundary mesh entities is imposed, regardless ofwhether it causes self-intersections. Then, the minimization proceeds in iterativelyupdating N towards lower energy configurations, that is: N k → N k +1 , such that E k → E k +1 < E k . The minimization is achieved via a gradient- and Hessian-basedmethod, and utilizes a relaxation strategy, whereby a local optimization problem issolved 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 posi-tion of the node under consideration. This requires the solution of a smaller and lesscomputationally expensive problems, compared to the global problem that would beotherwise necessary to solve. The minimization of equation (2.5) can then be recastinto the minimization of a set of energies, each of which belongs to a node and thesubset of elements the node influences. The stopping criterion used to find the globalminimum based on the set of local energies, is || N k +1 − N k || ∞ /L < (cid:15) , where we used (cid:15) = 10 − , and || s || ∞ = max {| s | , . . . , | s n |} is the uniform norm.This process also requires the computation of the gradient of E ( φ ), with respectto the movement of node positions. As a first approach, these gradients can beapproximated by use of finite difference stencils. However, our testing indicates thatthis approach is both expensive and sensitive to the size of the stencil. Instead, wemake use of an analytic expression for the desired gradients. Ultimately, applicationsof the chain rule therefore require the same derivatives for φ I and φ M . Although theideal 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 thederivatives can be obtained, which are given in the appendix of [80].Figure 4 shows the high-order surface mesh of the Rp1 car, resulting from ourmeshing process. It is possible to appreciate the high-order elements tessellating theCAD surface. We note that the spokes within the wheel shown in figure 4 werenot included in the simulations presented in section 4 to avoid the requirement for amoving mesh.As a final remark, we point out that, all those nodes on the boundary surface meshthat belonged to the “exclusion set” are left straight-sided. This leads to a small butnot insignificant proportion of the surface mesh entities that will be straight-sided asopposed to curved, and could have significant impact on the geometric accuracy of themesh and hence the solution. However, in the results that we present in section 4, weobserve little noticeable impact, since often these excluded elements were in locationswhere multiple CAD surfaces met but not typically in places of interest to the flowphysics. Obviously this cannot be always guaranteed to be true. A full investigationinto this has not been conducted yet, but this work proves that high-order meshes can2 G. MENGALDO, ET AL.
Fig. 4 . Details of the Rp1 high-order mesh produced using the high-order meshing workflowpresented in section 2. We note that the spokes in the wheel were not included in the final simulationsto reduce complexity. be produced on industrial models and sensible results obtained with little meshingeffort. Indeed the final projection meshing pipeline only required one execution of thelinear mesh generation and
NekMesh , with no need for repetition with either systemor the CAD.
3. Numerical methodology.
The model underlying the simulations carriedout in this study is constituted by the incompressible Navier-Stokes equations for aNewtonian fluid(3.1a) ∂ u ∂t + u · ∇ u = −∇ p + ν ∇ u + f (3.1b) ∇ · u = 0 , where u = ( u, v, w ) is the velocity vector, p denotes the pressure, ν is the kinematicviscosity, and f is a generic forcing term. This set of equations are particularly suitedto accurately model the flow around a road car, as compressibility effects are deemedto be unimportant at the speeds reached. For diffusers and other heated parts of thecar, the compressible Navier-Stokes equations are better suited. This was however outof the scope of the current investigation, that was mainly focussed on the externalaerodynamics of the car.The system of equations (3.1) and its numerical discretization has been imple-mented in the open-source spectral/ hp element library Nektar++ . The numerical dis-cretization uses an operator splitting scheme (sometimes referred to as a velocity cor-rection scheme) [36, 29] in combination with a consistent boundary condition for thepressure Poisson equation. The spatial discretisation follows a high-order continuousGalerkin formulation [38], where the elemental basis functions are typically modifiedJacobi or Lagrange polynomials defined on a set of Gauss-Lobatto-Legendre quad-rature points. The time-integration is instead achieved via a second-order implicit-explicit time-stepping strategy [82].In order to successfully achieve the challenging high-Reynolds number LES de-scribed in section 4, we needed to significantly enhance the numerics implemented in
NDUSTRY-RELEVANT IMPLICIT LARGE-EDDY SIMULATION Nektar++ , especially in terms of (i) numerical stability, (ii) regularization techniquesfor under-resolved scales, and (iii) computational efficiency. These three key buildingblocks are discussed in the following.
Numerical stability was a key aspect that we neededto carefully take into account in order to successfully simulate the road car presentedin section 4, at the Reynolds numbers required by industry. To achieve this, weadopted polynomial dealiasing via consistent integration [48], combined with spectralvanishing viscosity (SVV) [35, 40, 50]. The former helps suppressing aliasing-driveninstabilities that are due to inexact numerical integration of the equations, while thelatter targets all those flow scales that are necessarily under-resolved in high-Reynoldsnumber simulations. While SVV primary target was not the numerical stability of thesimulations, 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 insection 3.2, where we discuss its use to accurately treat under-resolved scales. In thissection, we focus on polynomial dealiasing that was a key component without whichit was otherwise not possible to obtain stable simulations.Indeed, the main issue when dealing with element-wise integrations typical ofspectral/ hp methods, is the inexact numerical integration of the discrete terms of theequations. In fact, it is well known that there is a minimum number of quadraturepoints Q min , that is required to produce exact results for integrals of polynomials ofany given degree. This minimum number of quadrature points depends on the type ofnumerical quadrature adopted and on the set of quadrature points used. In our case,we employed Gaussian quadrature on a set of Gauss-Lobatto-Legendre quadraturepoints [38]. This type of numerical integration guarantees that the minimum numberof 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 spaceof polynomials. It is then possible to calculate the number of quadrature pointsnecessary 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 thenonlinearity increases. Polynomial order P Q min [ u ( ξ )] ∈ P P Q ≥ P + 3 / u ( ξ )] ∈ P P Q ≥ P/ / u ( ξ )] ∈ P P Q ≥ P + 3 / Table 1
Number of GLL quadrature points for the GLL quadrature to be exact up to machine precisionas a function of the polynomial order of the integrand
If we now consider the incompressible Naver-Stokes equations (3.1), and we splitthem into linear and nonlinear terms as follows:(3.2) dudt = M − ( L + N ) , where M is the mass matrix that arises from the numerical discretization of the equa-tions, L = −∇ p + ν ∇ u + f contains the linear terms, and N = − u · ∇ u contains theconvective term that is a quadratic nonlinearity, we can identify the Q min requiredfor exact integration. Given that the spatial discretization we employ is a continuous4 G. MENGALDO, ET AL.
Galerkin method, we are interested in the integration of L inner products of twopolynomials ( u p , u q ). On the one hand, to calculate the linear terms, it is necessaryto calculate [ u ( ξ )] ∈ P P (i.e. the inner product of the solution with respect to thepolynomial basis). Therefore, for exact integration, we need Q min = P + 2 Gauss-Lobatto-Legendre quadrature points in one-dimension and a tensor product of pointsin higher dimensions. On the other hand, to calculate the nonlinear terms, we wish tocalculate [ u ( ξ )] ∈ P P , that corresponds to Q min = ( P +1) Gauss-Lobatto-Legendrequadrature 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 typicallypossible to exactly integrate the nonlinearities N present in the equations when deal-ing with mesh elements that are straight sided parallelograms, it is more complicatedwhen dealing with curved mesh elements. In fact, in the case of straight side paral-lelograms, the Jacobian of the mapping that defines the coordinates of the element isaffine with constant determinant. For curved elements however, this mapping is anisoparametric polynomial expansion, which when incorporated into integrands leadsto an additional source of aliasing error, called geometrical aliasing [48], that couldcouple with the underlying equation nonlinearities, thereby exacerbating numericalinstabilities.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 termof the numerical discretization (3.2) of (3.1), that is(3.3) dudt = M − (cid:12)(cid:12)(cid:12) Q (cid:18) L (cid:12)(cid:12)(cid:12) Q + N (cid:12)(cid:12)(cid:12) Q (cid:19) , The subscript Q in each term of equation (3.3) denotes a higher number of quadraturepoints than required for exact integration of quadratic linear terms, that is Q >P + 3 / Q > P/ /
2. We note that for simulations that are adequately resolved (e.g. lam-inar flows), aliasing usually does not affect numerical robustness. However, for under-resolved turbulence computations, this aliasing effect leads to a significant build-up oferror in high-frequency solution modes and usually causes the simulation to abruptlydiverge.There are additional tools to suppress aliasing-driven instabilities, that have notbeen explored in this work. These include the use of skew-symmetric forms of theconvective term of the underlying equations [2, 25, 83], 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 tosuppress aliasing-driven instabilities [35]. However, the diffusion required would besignificantly high, and could suppress important flow features, thereby underminingthe accuracy of the simulation. Additionally, because the SVV operator is relativelyisotropic, whereas dealiasing is anisotropic, the two do not completely overlap. Hence,
NDUSTRY-RELEVANT IMPLICIT LARGE-EDDY SIMULATION
15a reduced amount of regularisation can be achieved with dealiasing, thereby betterpreserving the solution accuracy.The use of the dealiasing techniques described in this section removes numericalintegration errors that could lead to numerical instabilities in under-resolved simula-tions; however, it does not address the under-resolved scales that are ubiquitous inhigh-Reynolds number simulations. These under-resolved scales, similarly to aliasingerrors, can be one of the main drivers of numerical instabilities and solution inac-curacies, and they exacerbate aliasing errors leading to the possible failure of thesimulation. The handling of under-resolved scales to prevent numerical instabilitiesand maintain solution accuracy is described next.
In order to treat under-resolved scales commonly present in high-Reynolds number simulations, such as theone 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 fromrecent dispersion/diffusion eigenanalysis [56, 49, 47, 50], and results in well-balanceddiffusion characteristics that allow for the accurate simulation of high Reynolds num-ber flow problems. Indeed, when stabilising a high-order formulation, such as thecontinuous Galerkin scheme adopted in this study, it is critical to achieve an adequatebalance between numerical diffusion and solution accuracy, especially for turbulentflow computations.The general idea behind SVV is that its diffusion can be designed to affect mainlythe element-wise polynomial modes of highest order within the spectral/ hp elementdiscretization. These polynomial modes comprise the most oscillatory part of the so-lution and represent the smallest turbulent scales captured by the degrees of freedomemployed. The SVV technique was first proposed by Tadmor [75] as a stabilisationstrategy for classical spectral methods in Fourier space. Essentially, an artificial dif-fusion term is added to the governing equations in the form(3.4) ν svv ∇ ( Q (cid:63) ∇ u ) = − ν svv (cid:88) κ κ ˆ Q κ ˆ u κ exp( iκx ) ,where ν svv is a constant representing SVV’s overall diffusion magnitude, κ is thewavenumber, ˆ u κ are the Fourier coefficients of the solution, and ˆ Q κ are the Fouriercoefficients of the SVV kernel that define how diffusion is to be distributed acrossFourier modes. We note that by setting ˆ Q κ ≡
1, the standard diffusion (Laplacian)operator is formally recovered. Traditional SVV kernel functions grow monotonicallyfrom zero, achieving ˆ Q κ = 1 only at the highest order mode. For that solutioncomponent, the physical kinematic viscosity ν effectively becomes ν + ν svv .Since Tadmor’s original work on spectral methods, SVV has been adapted tocontinuous spectral/ hp element methods and used in particular to stabilise under-resolved turbulence computations based on the continuous Galerkin method [35, 40,41, 43, 70]. In this approach the filtering properties of the SVV operator applied tothe continuous Galerkin (also referred to as CG hereafter) method are used in-lieu ofan explicit turbulence model, similarly to discontinuous Galerkin (also referred to asDG hereafter) methods, where the diffusion of the underlying numerical discretizationthat originates from the (upwind) Riemann fluxes act as a filter for the under-resolvedscales [12, 81, 53]. This setting, where artificial numerical features are used in place ofan explicit turbulence model, is commonly referred to as implicit LES, and it has beeninvestigated in the past few decades with various degrees of success [45, 55, 9, 10, 51].The key issues when modeling under-resolved scales using the diffusion properties6 G. MENGALDO, ET AL. of the underlying numerical discretization or through an SVV-like operator are lackof robustness and/or accuracy of the resulting implicit LES framework, and unclearconnections 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 [56, 49, 47, 50] to build anew SVV kernel. In particular, we used numerical dispersion and diffusion estimatesfor spectral/ hp element methods in wavenumber/frequency space provided by the so-called temporal ESA [56] to construct an SVV kernel Q that enforces SVV diffusionto be distributed monotonically across wavenumber space for the CG method. Wenote in fact that monotonic behavior in wavenumber space is not always guaranteeddue the lack of one-to-one correspondence between Fourier and polynomial modes .This resulted in what we called “power” CG-SVV kernel(3.5) ˆ Q ( p ) = (cid:18) pP (cid:19) P svv , 0 ≤ pP ≤ p is the polynomial mode index, P is the overall polynomial order, and P svv isa user-defined parameter that, when increased, further confines CG-SVV damping tothe smaller scales. This kernel also guarantees that raising P increases the resolutionpower per degree of freedom and is suitable for problems that assume periodic bound-ary conditions, as its construction is based on the temporal ESA framework [56]. It istherefore relevant to temporally evolving simulations versus spatially evolving ones,where inflow and outflow boundary conditions are typically used. While the “power”kernel defined in equation (3.5) solved a number of issues, still some potentially unde-sirable characteristics were present. These undesirable aspects were two-fold. On theone hand, non-smooth behavior of the diffusion and dispersion curves was observed foreither too small or too large values of ν svv . This could lead to unphysical representa-tion of the flow physics. On the other hand, the use of the temporal ESA frameworkdid not include spatially evolving flow problems (i.e. problem with inflow/outflowboundary conditions) that are commonly found in external aerodynamics. This couldtrigger spurious reflections of certain waves when they hit regions of mesh expan-sion (e.g. along turbulent wakes), eventually leading to numerical instabilities anddivergence of the simulation.In order to address the first aspect, we focused on the relevant non-dimensionalparameter to quantify SVV’s relative intensity is the so-called P´eclet number, namelyPe = vh/ν svv , that is based on measures of local velocity v and mesh spacing h .The adopted strategy was to maintain the P´eclet number fixed by making ν svv ∝ vh .This led to a CG-SVV scaling that is dimensionally equivalent to the upwind-baseddissipation of discontinuous spectral/ hp element methods.To address the second aspect, we pursued a complementary study [50] focused onthe so-called spatial ESA, which assumes inflow-/outflow-type boundary conditionsand is therefore more relevant to spatially developing problems, as the one encounteredin this paper.
Spatial
ESA had initially been conducted for discontinuous spectral/ hp element methods [49, 47] and proved insightful in explaining how spurious reflectionsoccurred in regions of mesh expansion depending on the Riemann flux employed.Note that for an incoming wave of given frequency ω , pure advection with velocity a should yield a single real-valued wavenumber κ = ω/a . However, in a numericalsetting, two complex-valued wavenumbers, κ p (physical) and κ u (unphysical), appear Note that common practice was to directly apply Fourier-type kernel entries ˆ Q κ upon theorthogonal hierarchical set of element-wise polynomial modes.NDUSTRY-RELEVANT IMPLICIT LARGE-EDDY SIMULATION (cid:61) ( κ p ) > u ∝ exp[ i ( κx − ωt )] = exp[ −(cid:61) ( κ ) x ] exp { i [ (cid:60) ( κ ) x − ωt ] } ,where (cid:61) and (cid:60) correspond to the imaginary and real part, respectively. For theunphysical reflected mode, it is instead necessary that (cid:61) ( κ u ) is negative, indicatingdamping during backward propagation (∆ x < I ( κ u ). This allows for spurious waves to survive for along time while interacting with incoming turbulence, eventually causing numericalinstability. Hence, a new optimisation procedure was performed for the CG-SVVoperator in order to match the diffusion levels of CG to those of standard upwind DGat appropriate polynomial orders, this time with a constraint to enforce a sufficientlystrong diffusion for the reflected eigenmodes. DG’s numerical diffusion was taken as areference because DG is regarded as having an adequate balance between dissipationand accuracy. The resulting dispersion and diffusion plots are shown in figure 5 forvarious orders, with the curves for the spurious modes at the bottom part of thegraphs. Once again, the newly designed kernel for the CG-SVV operator, that wenamed “DG” kernel given its connections to the DG method in terms of diffusion anddispersion characteristics, maintains the P´eclet number fixed and guarantees thatraising P increases the resolution power per degree of freedom. Fig. 5 . Numerical dispersion (left) and diffusion (right) plots for CG with the “DG kernel”SVV for various polynomial orders. The unphysical eigenmodes appear on the bottom curves of thegraphs. Figure taken from Moura et al. [50].
In order to validate the novel CG-SVV operator based on the “DG” kernel, weused a simple two-dimensional test case that was designed to mimic spatially devel-oping grid turbulence [49, 50]. It consists of a rectangular domain with prescribedunsteady flow at the inlet and symmetry (free-slip) conditions on the sides. Thisdomain was meshed with quadrilaterals in a structured grid pattern featuring tworegions. The upstream region had square-shaped elements, whereas the downstream8
G. MENGALDO, ET AL. region had a four times larger streamwise mesh spacing. The interface between thesetwo regions marked a potentially unstable threshold zone across which mesh spacingchanged abruptly. Results obtained without SVV for P = 8 and n = 11 elementsin the crossflow direction are shown in figure 6. In particular, the vorticity contourplot obtained without SVV, upon close inspection, shows the presence of spuriousreflections 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, higherReynolds numbers were required for the simulation to diverge, while for higher poly-nomial orders, the test cases started to diverge at correspondingly lower Reynoldsnumbers. When the novel DG-based CG-SVV operator was added, however, up topolynomial order P = 8, all cases remained stable no matter how high the Reynoldsnumber tested. Furthermore, no spurious reflections have been noticed. And althoughsolutions at lower orders ( P = 3 in particular) might have seemed excessively regu-larised by the “DG” kernel SVV, superior resolution power was effectively recoveredat higher orders, such as P = 4 and above [52]. Fig. 6 . Vorticity contours of case P = 8 , n = 11 , at Re = 1000 with (bottom plot) andwithout (top plot) the novel CG-SVV operator in presence of streamwise mesh coarsening. Closeexamination shows how the novel CG-SVV kernel is able to suppress spurious reflections originatingwhere the mesh changes resolution. Figure readapted from Moura et al. [50]. The current best practice in
N ektar has been to use the “DG” kernel in generalwhile employing the “power” kernel only along homogeneous directions in case peri-odic boundary conditions are present. In either case, SVV’s base coefficient takes thefinal form(3.7) ν svv = (Pe ∗ ) − v h/P ,where Pe ∗ is a reference constant Pecl´et number (typically of order one), v and h arerespectively local measures of advection velocity and mesh size (along the relevantdirections of flow propagation), and P is the chosen polynomial order. The entriesdefining the “power” kernel are those given in (3.5), whereas those of the “DG” kernelhave no explicit formula (as the output of an optimisation procedure), having theirvalues given in table form for each order P in [50]. Furthermore, we note that SVV’skernel operation from (3.4) is accommodated in polynomial space as(3.8) Q (cid:63) ∇ u ≈ R − Q R ˆg ,in which ˆg denotes the vector of element-wise coefficients in the polynomial expansionof ∇ u , while R is a “rotation” matrix that projects ˆg onto an orthogonal (Legendre) NDUSTRY-RELEVANT IMPLICIT LARGE-EDDY SIMULATION Q is the diagonal matrix whose entries are given by the chosenSVV kernel function. Therefore, the proposed CG-SVV operator has the final form(3.9) ν svv ∇ ( Q (cid:63) ∇ u ) ≈ (Pe ∗ ) − v hP ∇ (cid:0) R − Q R ˆg (cid:1) .This CG-SVV operator has been implemented in
Nektar++ and extensively tested indifferent test cases of increasing complexity, always displaying strong robustness andpreventing numerical instabilities. This approach offers a decent balance between ro-bustness and accuracy for high-order discretizations, being effectively a parameter-freestrategy from the practitioners’s perspective and allowing for under-resolved turbu-lence simulations over complex geometries at high Reynolds numbers. This approach,in conjunction with the dealiasing techniques described in section 3.1, was used toobtain the numerical simulations depicted in section 4. Without the concurrent use ofdealiasing and SVV regularization, the numerical simulations presented in this paperwould have diverged due to numerical instabilities.
One of the advantages of higher order poly-nomial expansions is their ability to exploit high arithmetic intensities; that is, agreater number of floating point operations is performed for each byte of data re-trieved from memory at higher orders than at lower orders. Although this naturallymakes operator evaluations at higher order more expensive per degree of freedom, onmodern computing hardware this effect frequently tends to be masked, as the costof transferring data from memory to the compute unit tends to be far greater thanthat of performing floating point operations. In essence, this can make floating pointoperations nearly ‘free’ compared to the transfer of data.However, the efficient implementation of the building-block finite element opera-tors 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 requirescareful consideration. Even small changes in polynomial order can yield a large ef-fect on performance due to the many dependent factors, including both the choice ofhardware and the problem under consideration. Our recent work in this area [57] hasfocused on the development of techniques that can automatically determine the mostefficient kernels for these operators at runtime. These kernels are aimed to targetboth clear performance factors (such as operator type and polynomial order), as wellas more subtle aspects, such as whether the choice of basis and whether mathemat-ical techniques such as sum-factorisation can be leveraged to increase performance,and the grouping of elements with equivalent numerical properties to improve per-formance evaluation. The use of these kernels therefore allows an efficient on-nodeimplementation across a broad range of polynomial orders.For parallelisation between nodes,
Nektar++ uses a parallelisation model based ondistributed MPI. The METIS library [39] is first used to partition the domain basedon the dual graph representation of the mesh. Since, for industrial geometries, themesh generation process yields a hybrid mesh of prismatic and tetrahedral elements,these are weighted appropriately in the partitioning process, so as to approximatelyequally distribute the work between processors. Both the pressure Poisson equationand the three scalar velocity correction Helmholtz equations that comprise the oper-ator splitting scheme are solved implicitly, leading to the need for efficient solvers forlarge, distributed linear algebraic positive definite systems. To achieve this, we applyan iterative preconditioned conjugate gradient (PCG) method for each scalar system.We note, however, that the usual approach taken at lower orders of representing this0
G. MENGALDO, ET AL. 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 itera-tions, we implement this as many smaller, locally dense elemental matrices, combinedwith a sparse gather operation to apply C connectivity. The gather operation usesthe efficient Gslib library in parallel, based on the code Nek5000 [18]. Our testshave demonstrated that this combination allows both performance and scalability ofthe PCG solver up to 10 cores at high levels of efficiency. Finally, we discuss pre-conditioning of these methods. All four PCG solvers use an effective and scalablelow-energy basis preconditioner [72, 28] and, optionally, projected initial conditionsbased on storage from solutions at previous timesteps [20]. The solution of the pres-sure Poisson equation, however, is particularly ill-conditioned, and therefore requiresthe use of coarse-space preconditioning, which requires the inversion of the linear fi-nite element representation of the mesh. Currently, the coarse space preconditioneruses a direct parallel inverse method developed by Fischer and Tufo [78].The combination of all of these intra- and inter-node performance optimisations,together with a robust meshing pipeline and stabilisation strategies, has yielded an ef-ficient solver at high orders that is capable of tackling complex configurations requiredfor industrially relevant problems at superior computational performance.
4. High-order iLES for the aerodynamic (re)design of the ElementalRp1 road car.
The developments in terms of high-order mesh generation and nu-merics implemented in
NekMesh and
Nektar++ that were described in sections 2 and3 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 Ele-mental Cars since June 2016, is a high-performance, light-weight, road-legal vehiclethat weighs 580 kg. It is equipped with a 320 bhp engine and accelerates from 0to 60 mph in just 2.6 seconds. Its current configuration, depicted in Figure 7 alongwith its digital twin, delivers 400kg of downforce at 150 mph. The main targets ofthe high-fidelity LES simulations performed were to analyze and improve the externalaerodynamic 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 isthe baseline design, while D2 and D3 are two design upgrades that were implementedto improve aerodynamic performance and pilot comfort. This study was a collab-orative 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 theopen-source library
Nektar++ . In the following, we describe the configuration of thesimulations and we discuss qualitatively the results obtained.(a) (b)
Fig. 7 . The original design of the Elemental Rp1 road car, a high performance, street legaltrack car which delivered 400 kg of downforce at 150 mph.
NDUSTRY-RELEVANT IMPLICIT LARGE-EDDY SIMULATION Thissection outlines the details of the meshes adopted along with the strategies used togenerate them (section 4.1.1), and discusses the setup of the simulations performed(section 4.1.2).
The hugely complicated nature of the geometry un-der consideration, which involves hundreds of CAD surfaces and thousands of edges,presented a significant test for the mesh generation process outlined in section 2. Thehigh-order meshes used in this study ranged from approximately 2.2 million to 4 mil-lion elements, and they are comprised of a combination of triangular prisms in thenear-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 methoddescribed in section 2.3, for the surface mesh generation stage. The first analyticcurving method was applied to the first design, namely D1, where the creation ofa viable mesh took approximately 4 months of iterations between CAD healing andhigh-order meshing. The final mesh was composed of 2.2 million elements, including600000 prisms within the boundary-layer mesh. The second analytic curving methodwas instead applied to the second design, namely D2, that was constructed basedon the results from D1. The mesh for D2 contained 2.6 million elements, including700000 boundary layer prisms. While, both analytic curving strategies were effectivein generating viable meshes for D1 and D2, they were deemed excessively user-drivenand time-consuming. The mesh for D3 was obtained with what we named the projec-tion curving method, described in section 2.3. Indeed, this method, along with theuse of
CADfix for geometry repair,
Star-CCM+ for linear and boundary-layer meshgeneration and
NekMesh for high-order mesh generation was the least time-consumingand user-driven. It allowed the mesh to be created almost without user effort, andit required one single execution of linear mesh generation via
Star-CCM+ and onesingle execution of
NekMesh . No repetitive cycles of CAD healing were required. Theoverall meshing process for D3 produced a mesh with approximately 3.5 million el-ements, including 1.24 million boundary-layer prisms, and it took approximately 1week. Although the meshing strategy for D3 is the most robust and efficient, weshould point out that allows for inaccuracies in the geometry description, as a smallpercentage of surface elements are included in the “exclusion set”, thereby remainingstraight-sided. However, in the results that we present in section 4.2.2, no obviousunphysical results could be noticed, and the flow physics in those regions of the carthat were geometrically identical between D2 and D3 was indistinguishable.Figure 8 shows the three surface meshes for D1, D2, and D3 obtained via thethree different strategies just outlined. It is possible to note that the surface meshfor D2 and D3 has a smoother element distribution than D1, due to the enhancedtriangulation strategies adopted, namely the second analytic curving , and projectioncurving . We remark that for all the simulations carried out we did not include thespokes that are instead shown in figure 4.
The simulation setup required to analyze the aero-dynamic of the Rp1 car was challenging due to the high Reynolds number, Re = 10 ,based on the length of the car L , and to the complex geometries involved. In orderto obtain stable yet accurate simulations, we adopted some simplifications that arecommonly adopted when performing the aerodynamic car simulations.The Rp1’s unique double diffuser design and the the road-tyre interaction were2 G. MENGALDO, ET AL. (a)
High-order surface mesh (b)
Top down view (c)
Side on view
Figure 6.2. (a) shows the high-order surface mesh, (b) and (c) two images on Cp = 0isocontours coloured in pressure, one with a top down view the other sideon. In each case the half car on the on the left is the design 1 geometryand corresponding flow and design 2 on the right. redesigned console area in front of the driver. The D2 isocontour shows clearly thatthese structures now pass cleanly over the drivers head and are in somewhat lessnoisy in terms of the substructures. Secondly, the roll hoop produced significantlymore drag and separation than predicted. D2 had a redesigned roll hoop with an Design 1 (D1) Design 2 (D2) (a)
High-order surface mesh (b)
Top down view (c)
Side on view
Figure 6.5. (a) shows the high-order surface mesh, (b) and (c) two images on Cp = 0isocontours coloured in pressure, one with a top down view the other sideon. In each case the half car on the on the left is the design 1 geometryand corresponding flow and design 3 on the right. alteration in ride height. The results shown in figure 6.5, are consistent with the D2simulation and are shown alongside the D1 results. Again the trends in increasing Design 1 (D1) Design 3 (D3) (a)
High-order surface mesh (b)
Top down view (c)
Side on view
Figure 6.2. (a) shows the high-order surface mesh, (b) and (c) two images on Cp = 0isocontours coloured in pressure, one with a top down view the other sideon. In each case the half car on the on the left is the design 1 geometryand corresponding flow and design 2 on the right. redesigned console area in front of the driver. The D2 isocontour shows clearly thatthese structures now pass cleanly over the drivers head and are in somewhat lessnoisy in terms of the substructures. Secondly, the roll hoop produced significantlymore drag and separation than predicted. D2 had a redesigned roll hoop with an Design 1 (D1) Design 2 (D2) (a)
High-order surface mesh (b)
Top down view (c)
Side on view
Figure 6.5. (a) shows the high-order surface mesh, (b) and (c) two images on Cp = 0isocontours coloured in pressure, one with a top down view the other sideon. In each case the half car on the on the left is the design 1 geometryand corresponding flow and design 3 on the right. alteration in ride height. The results shown in figure 6.5, are consistent with the D2simulation and are shown alongside the D1 results. Again the trends in increasing Design 1 (D1) Design 3 (D3) (a) (b)
Fig. 8 . Surface meshes of the three car designs considered, D1, D2, and D3. Figure (a) showsa comparison between D1 and D2, while figure (b) shows a comparison between D1 and D3. It ispossible to notice a smoother distribution of surface elements for both D2 and D3 when comparedto D1. This is due to the different approaches adopted for the surface mesh generation, where forD1 we used the first analytic curving, while for D2 and D3 we used the second analytic curving andprojection curving, respectively. carefully modeled with some simplifications that allowed a significantly easier simula-tion setting while not undermining the accuracy of the results. Rp1’s double diffuserdesign has intakes and outlets on the body of the car. To model these accurately, weextracted the corresponding surfaces and constructed an appropriate velocity bound-ary field, by enforcing a prescribed average normal velocity on each surface. Theboundary condition has to also satisfy the no-slip condition at the edges of the inletsurfaces and this was imposed by solving a Helmholtz problem where the positiveHelmholtz 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]. Inaddition, the road-tyre interaction was modeled using a contact patch along with ro-tational boundary conditions to a simplified tyre geometry. This allowed efficientlymodeling the rotation of the tyre without the need for a moving mesh. Finally, wesimulated only half of the car, by using a symmetry plane. These simplificationsallowed a faster simulation setup and simulation runtime, without altering the flowfeatures 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 . The initial conditions wereimposed impulsively so that u/U = 1, where U is a reference velocity, and all non-moving surfaces are treated as no-slip boundaries. At the outflow, we imposed theboundary condition developed by Dong et al. [14], which balances the kinetic energyinflux through the outflow boundary condition to prevent instability. The Reynoldsnumber is then gradually increased by a factor of ten until the target Reynolds number, Re = 10 , is reached. Between increases, the flow is evolved for two convective timeunits t c = L/U to allow the flow to adapt to the new flow physics dictated by thehigher Reynolds number. We note that the Reynolds number used in this study, Re = 10 , is still lower than the experimental value of approximately Re = 2 × . NDUSTRY-RELEVANT IMPLICIT LARGE-EDDY SIMULATION hp elementmethods.The simulations were run at a polynomial order 4 (fifth order of accuracy) witha second order implicit-explicit time-stepping. To increase the accuracy of the inte-gration, a seventh order of quadrature was used, that consisted of the so-called globalconsistent integration introduced in Mengaldo et al. [48], and described in section 3.1.This allowed to mitigate aliasing-driven instabilities. We note that the use of a higherquadrature meant that the meshes were generated at a higher order of P = 7 to ensurethe elements were valid with this quadrature. In addition, the SVV kernel presentedin section 3.2 was used to filter the under-resolved scales and enhance the stabilityof the simulations. Without the use of both, consistent integration and SVV was notpossible to obtain numerically stable calculations.In terms of computational requirements, the simulation were performed in parallelusing 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. Thehigh-fidelity flow simulations performed allowed unique insights into the aerodynamicperformance of the Rp1’s complex geometry, highlighting a number of importantregions for aerodynamic optimisation. These results are discussed next. In this section, we discuss the results obtained on the three cardesigns introduced earlier, D1, D2 and D3. In particular, in section 4.2.1, we presenta comparison between D1 and D2, while in section 4.2.2, we present a comparisonbetween D1 and D3. The results focus on iso-contour quantities that allowed toidentify important flow features suitable to optimization.
The first simulation performed was for design 1(D1), that is the first baseline model of the car that was produced, followed by thesimulation of design 2 (D2). The latter design, D2, implements an upgrade of theaerodynamic package with respect to D1. Figure 9 shows the meshes and associatedresults 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 flowresults. 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-contoursof the total pressure coefficient Cp = 1 − ( u/u ∞ ) = 0, and they are colored aspressure. From the high-order D1 simulation two key findings were made, neither ofwhich were identified in low-order RANS simulations [76]. Firstly, there appeared tobe 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 thewindscreen, and the issue was fixed in D2 with a redesigned console area in front ofthe driver, significantly enhanced the driving experience. The D2 iso-contour showsclearly that these structures now pass cleanly over the drivers head and are less noisyin terms of small-scale vortices (figures 9-(b) and 9-(c), for D2). Secondly, the rollhoop produced significantly more drag and separation than predicted. D2 had aredesigned roll hoop with an aerofoil profile as opposed to a circular cylinder. Thisaerofoil profile was slightly angled to help control the flow over the new Gurney flapat 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 dragbetween the D1 and D2 simulations was well predicted by the high-order simulationsand matched the trend in the RANS results [76].4 G. MENGALDO, ET AL.
Design 1 (D1) Design 2 (D2) Design 1 (D1) Design 2 (D2)
Design 1 (D1) Design 2 (D2) Design 1 (D1) Design 2 (D2) (a) (b)
Design 1 (D1) Design 2 (D2)
Pressure0.0-1.0-0.5-0.25-0.75
Wake from windscreen on D1 (c)
Fig. 9 . Comparison between D1 and D2, in terms of mesh [figure (a)], and flow physics [figures-(b) and -(c)]. The images depicting the flow physics shows the total pressure coefficient for Cp = 0 and the iso-contours are colored in pressure. The final simulation performed in this study is fordesign 3 (D3), which is a full aerodynamic upgrade over the previous two designs, D1and D2. This car is designed to achieve extremely high levels of downforce specificallyfor track racing. This design includes a fully redesigned floor, front splitter and theaddition of a rear wing. The ride height has also been altered, raising the car at therear for increased diffuser performance and lowering the front of the car to increasein-ground effect of the front splitter.Figure 10, similarly to figure 9, shows again the total pressure coefficient, coloredby pressure. The results are consistent with the D2 simulation and are shown alongsidethe D1 results. The trends in increasing downforce were well predicted by the high-order simulation. In figure 10, we note that it is possible to observe a noticeable offsetin 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 preservingmesh validity versus geometrical accuracy, thereby leaving a non-negligible numberof elements straight-sided as opposed to curvilinear, the flow physics is in agreementwith D1 and D2, and no unphysical results could be readily noticed. However a morecomprehensive study should be undertaken in order to draw stronger conclusions on
NDUSTRY-RELEVANT IMPLICIT LARGE-EDDY SIMULATION Design 1 (D1) Design 3 (D3) Design 1 (D1) Design 3 (D3)
Design 1 (D1) Design 3 (D3) Design 1 (D1) Design 3 (D3) (a) (b)
Design 1 (D1) Design 3 (D3)
Pressure0.0-1.0-0.5-0.25-0.75 (c)
Fig. 10 . Comparison between D1 and D3, in terms of mesh [figure (a)], and flow physics[figures (b) and (c)]. The images depicting the flow physics shows the total pressure coefficient for Cp = 0 and the iso-contours are colored in pressure. whether the compromised geometric accuracy of our most robust meshing pipeline isa 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, andsome 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 meshgeneration 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 thatconstitute the geometry of the car, and to avoid using an excessive number of linearelements in proximity of curved geometrical features. Indeed, the meshing frameworkhas proven reliable and viable for industrial geometries, and it was a key enabler tosuccessfully accomplish the simulations presented in section 4.Along with the high-order meshing framework, the novel numerical technologiesdescribed in this section 3 were essential to successfully produce the CFD workflowdescribed in this paper. More specifically, the consistent integration for mitigatingaliasing-driven instabilities in conjunction with the CG-SVV based on the “DG” ker-6
G. MENGALDO, ET AL. nel for treating under-resolved scales were the keys to obtain stable simulations andaccurate simulations. The use of the novel CG-SVV approach, in particular, wasinstrumental to produce high-fidelity results. In fact, the CG-SVV operator is basedon ESA’s estimates for spectral/ hp element methods, where numerical dispersion anddiffusion characteristics are quantified in wavenumber/frequency space, as in clas-sical von Neumann analyses. Although mostly feasible for simple model problems,such as one-dimensional linear advection, ESA’s estimates for spectral/ hp methodshave been shown to hold more generally for nonlinear problems and even turbulencesimulations [54, 52, 17]. These estimates quantify how the numerical error (diffu-sion in particular) increases at large wavenumbers/frequencies, affecting especiallythe smallest, poorly- and under-resolved scales. When interpreted alongside implicitLES experiments, ESA reveals that spectral/ hp discretisations of higher polynomialorders perform best because they affect the large and intermediate flow scales as littleas possible, while gradually introducing a sufficiently strong diffusion to regularisethe smaller scales. This way, the governing equations and thus the turbulence physicsare 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 onhow numerical diffusion acts as a turbulence model [45, 55, 9, 10, 51]. However theESA framework was deemed the most promising route to achieve high-fidelity implicitLES results.The simulations in section 4 represent the first fifth-order accurate results of a realautomotive geometry by means of implicit LES technologies using spectral/ hp elementmethods. They represent one of the first attempts at bridging the gap between appliedmathematics academic research and industrial applications within the context of CFD. Acknowledgments.
This work has benefited from the technical support anddesign insight from Elemental Cars and from access to large HPC resources providedby Hewlett Packard Enterprise. We also would like to thank Mark Gammon andhis group at ITI Global for providing us with a license to the program CADfix andits CFI interface, and for their technical help and support. Partial financial supportwas provided by EPSRC under the Platform Grant PRISM: Platform for Research InSimulation Methods (EP/R029423/1) and by the EU Horizon 2020 project ExaFLOW(grant 671571). M. Turner acknowledges the support of EPSRC and Airbus underan Industrial CASE studentship. We also acknowledge support from Imperial CollegeHigh Performance Computing Service.
REFERENCES[1]
J. L. Bentley , Multidimensional binary search trees used for associative searching , Commu-nications of the ACM, 18 (1975), pp. 509–517.[2]
G. A. Blaisdell, E. T. Spyropoulos, and J. H. Qin , The effect of the formulation of non-linear terms on aliasing errors in spectral methods , Applied Numerical Mathematics, 21(1996), pp. 207–219.[3]
J. H. Bucklow, R. Fairey, and M. R. Gammon , An automated workflow for high qualitycfd meshing using the 3d medial object , in 23rd AIAA Computational Fluid DynamicsConference, 2017, p. 3454.[4]
R. H. Bush, T. S. Chyczewski, K. Duraisamy, B. Eisfeld, C. L. Rumsey, and B. R. Smith , Recommendations for future efforts in RANS modeling and simulation , in AIAA Scitech2019 Forum, 2019, p. 0317.[5]
C. D. Cantwell, D. Moxey, A. Comerford, A. Bolis, G. Rocco, G. Mengaldo, D. DeGrazia, S. Yakovlev, J. Lombard, D. Ekelschot, B. Jordi, H. Xu, Y. Mohamied,C. Eskilsson, B. Nelson, P. Vos, C. Biotto, R. M. Kirby, and S. J. Sherwin , Nek-tar++: An open-source spectral/ hp element framework , Computer Physics Communica-NDUSTRY-RELEVANT IMPLICIT LARGE-EDDY SIMULATION tions, 192 (2015), pp. 205–219.[6] C. D. Cantwell, S. Yakovlev, R. M. Kirby, N. S. Peters, and S. J. Sherwin , High-orderspectral/ hp element discretisation for reaction–diffusion problems on surfaces: Applicationto cardiac electrophysiology , Journal of computational physics, 257 (2014), pp. 813–829.[7] G. Castiglioni and J. A. Domaradzki , A numerical dissipation rate and viscosity in flowsimulations with realistic geometry using low-order compressible Navier–Stokes solvers ,Computers & Fluids, 119 (2015), pp. 37–46.[8]
Centaur
J. B. Chapelier, G. Lodato, and A. Jameson , A study on the numerical dissipation of thespectral difference method for freely decaying and wall-bounded turbulence , Computers &Fluids, 139 (2016), pp. 261–280.[10]
T. Dairay, E. Lamballais, S. Laizet, and J. C. Vassilicos , Numerical dissipation vs subgrid-scale modelling for large eddy simulation , Journal of Computational Physics, 337 (2017),pp. 252–274.[11]
F. Dassi, A. Mola, and H. Si , Curvature-adapted remeshing of CAD surfaces , Procedia En-gineering, 82 (2014), pp. 253–265.[12]
C. C. de Wiart, K. Hillewaert, L. Bricteux, and G. Winckelmans , Implicit LES of freeand wall-bounded turbulent flows based on the discontinuous Galerkin/symmetric interiorpenalty method , International Journal for Numerical Methods in Fluids, 78 (2015), pp. 335–354.[13]
S. Dey, R. M. O’Bara, and M. S. Shephard , Towards curvilinear meshing in 3D: the caseof quadratic simplices , Computer-Aided Design, 33 (2001), pp. 199–209.[14]
S. Dong, G. E. Karniadakis, and C. Chryssostomidis , A robust and accurate outflowboundary condition for incompressible flow simulations on severely-truncated unboundeddomains , Journal of Computational Physics, 261 (2014), pp. 83–105, https://doi.org/10.1016/j.jcp.2013.12.042.[15]
Elemental Cars , Elemental cars . http://elementalcars.co.uk/, 2018.[16]
Elemental Cars , The Elemental RP1 . http://elementalcars.co.uk/the-rp1/, 2018.[17]
P. Fernandez, R. C. Moura, G. Mengaldo, and J. Peraire , Non-modal analysis of spectralelement methods: towards accurate and robust large-eddy simulations , Computer Methodsin Applied Mechanics and Engineering, 346 (2019), pp. 43–62.[18]
P. Fischer, J. Kruse, J. Mullen, H. Tufo, J. Lottes, and S. Kerkemeier , Nek5000: Opensource spectral element CFD solver . https://nek5000.mcs.anl.gov/index.php/MainPage,2008. Argonne National Laboratory, Mathematics and Computer Science Division, Ar-gonne, IL.[19]
P. Fischer, J. Mullen, et al. , Filter-based stabilization of spectral element methods , ComptesRendus de l’Academie des Sciences Series I Mathematics, 332 (2001), pp. 265–270.[20]
P. F. Fischer , Projection techniques for iterative solution of Ax = b with successive right-handsides , Computer Methods in Applied Mechanics and Engineering, 163 (1998), pp. 193–204.[21] P. F. Fischer, G. W. Kruse, and F. Loth , Spectral element methods for transitional flowsin complex geometries , Journal of Scientific Computing, 17 (2002), pp. 81–98.[22]
M. Fortunato and P.-O. Persson , High-order unstructured curved mesh generation usingthe Winslow equations , Journal of Computational Physics, 307 (2016), pp. 1–14.[23]
C. Fu, M. Uddin, and A. C. Robinson , Turbulence modeling effects on the CFD predic-tions of flow over a NASCAR Gen 6 racecar , Journal of Wind Engineering and IndustrialAerodynamics, 176 (2018), pp. 98–111.[24]
A. Gargallo-Peir´o, X. Roca, J. Peraire, and J. Sarrate , Optimization of a regularizeddistortion measure to generate curved high-order unstructured tetrahedral meshes , Inter-national Journal for Numerical Methods in Engineering, 103 (2015), pp. 342–363.[25]
G. J. Gassner , A skew-symmetric discontinuous Galerkin spectral element discretization andits relation to SBP-SAT finite difference methods , SIAM Journal on Scientific Computing,35 (2013), pp. A1233–A1253.[26]
G. J. Gassner and A. D. Beck , On the accuracy of high-order discretizations for underre-solved turbulence simulations , Theoretical and Computational Fluid Dynamics, 27 (2013),pp. 221–237.[27]
D. Gottlieb and J. S. Hesthaven , Spectral methods for hyperbolic problems , Journal ofComputational and Applied Mathematics, 128 (2001), pp. 83–131.[28]
L. Grinberg, D. Pekurovsky, S. J. Sherwin, and G. E. Karniadakis , Parallel performanceof the coarse space linear vertex solver and low energy basis preconditioner for spectral/ hp elements , Parallel Computing, 35 (2009), pp. 284–304.[29] J.-L. Guermond and J. Shen , Velocity-correction projection methods for incompressible flows ,SIAM Journal on Numerical Analysis, 41 (2003), pp. 112–134. G. MENGALDO, ET AL.[30]
R. Hartmann and T. Leicht , Generation of unstructured curvilinear grids and high-orderdiscontinuous Galerkin discretization applied to a 3D high-lift configuration , InternationalJournal for Numerical Methods in Fluids, 82 (2016), pp. 316–333.[31]
J. Hesthaven and R. Kirby , Filtering in Legendre spectral methods , Mathematics of Compu-tation, 77 (2008), pp. 1425–1452.[32]
F. Hindenlang, G. Gassner, T. Bolemann, and C. D. Munz , Unstructured high order gridsand their application in discontinuous Galerkin methods , in Conference Proceedings, VEuropean Conference on Computational Fluid Dynamics ECCOMAS CFD, 2010, pp. 1–8.[33]
T. J. R. Hughes, G. R. Feij´oo, L. Mazzei, and J.-B. Quincy , The variational multiscalemethoda paradigm for computational mechanics , Computer methods in applied mechanicsand engineering, 166 (1998), pp. 3–24.[34]
ITI-Global , CADfix: CAD translation, healing, repair, and transformation
G. S. Karamanos and G. E. Karniadakis , A spectral vanishing viscosity method for large-eddy simulations , Journal of Computational Physics, 163 (2000), pp. 22–50.[36]
G. E. Karniadakis, M. Israeli, and S. A. Orszag , High-order splitting methods for theincompressible Navier-Stokes equations , Journal of Computational Physics, 97 (1991),pp. 414–443.[37]
G. E. Karniadakis and S. A. Orszag , Nodes, modes and flow codes , Physics Today, 46 (1993),pp. 34–42.[38]
G. E. Karniadakis and S. J. Sherwin , Spectral/hp element methods for computational fluiddynamics , Oxford University Press, 2nd ed., 2005.[39]
G. Karypis and V. Kumar , A fast and high quality multilevel scheme for partitioning irregulargraphs , SIAM Journal on Scientific Computing, 2 (1998), pp. 359–392.[40]
R. M. Kirby and S. J. Sherwin , Stabilisation of spectral/hp element methods through spec-tral vanishing viscosity: application to fluid mechanics modelling , Computer Methods inApplied Mechanics and Engineering, 195 (2006), pp. 3128–3144.[41]
K. Koal, J. Stiller, and H. M. Blackburn , Adapting the spectral vanishing viscosity methodfor large-eddy simulations in cylindrical configurations , Journal of Computational Physics,231 (2012), pp. 3389–3405.[42]
LCS Fast , London computational solutions
J.-E. W. Lombard, D. Moxey, S. J. Sherwin, J. F. A. Hoessler, S. Dhandapani, and M. J.Taylor , Implicit large–eddy simulation of a wingtip vortex , AIAA Journal, 54 (2016),pp. 506–518.[44]
X.-J. Luo, M. S. Shephard, R. M. O’Bara, R. Nastasia, and M. W. Beall , Automaticp-version mesh generation for curved domains , Engineering with Computers, 20 (2004),pp. 273–285.[45]
L. G. Margolin and W. J. Rider , A rationale for implicit turbulence modelling , InternationalJournal for Numerical Methods in Fluids, 39 (2002), pp. 821–841.[46]
S. Marras, J. F. Kelly, F. X. Giraldo, and M. V´azquez , Variational multiscale stabi-lization of high-order spectral elements for the advection–diffusion equation , Journal ofComputational Physics, 231 (2012), pp. 7187–7213.[47]
G. Mengaldo, D. De Grazia, R. C. Moura, and S. J. Sherwin , Spatial eigensolution analy-sis of energy-stable flux reconstruction schemes and influence of the numerical flux onaccuracy and robustness , Journal of Computational Physics, 358 (2018), pp. 1–20.[48]
G. Mengaldo, D. De Grazia, D. Moxey, P. E. Vincent, and S. J. Sherwin , Dealiasingtechniques for high-order spectral element methods on regular and irregular grids , Journalof Computational Physics, 299 (2015), pp. 56–81.[49]
G. Mengaldo, R. C. Moura, B. Giralda, J. Peir´o, and S. J. Sherwin , Spatial eigensolu-tion analysis of discontinuous Galerkin schemes with practical insights for under-resolvedcomputations and implicit LES , Computers & Fluids, 169 (2018), pp. 349–364.[50]
R. Moura, M. Aman, J. Peir´o, and S. J. Sherwin , Spatial eigenanalysis of spectral/ hp continuous Galerkin schemes and their stabilisation via DG-mimicking spectral vanishingviscosity for high Reynolds number flows , Journal of Computational Physics, 406 (2020),p. 109112.[51] R. C. Moura, G. Mengaldo, J. Peir´o, and S. J. Sherwin , An LES setting for DG-based im-plicit LES with insights on dissipation and robustness , in Spectral and High Order Methodsfor Partial Differential Equations - ICOSAHOM 2016, M. Bittencourt, N. Dumont, andJ. S. Hesthaven, eds., Springer, 2017, pp. 161–173.[52]
R. C. Moura, G. Mengaldo, J. Peir´o, and S. J. Sherwin , On the eddy-resolving capabilityof high-order discontinuous Galerkin approaches to implicit LES / under-resolved DNS ofEuler turbulence , Journal of Computational Physics, 330 (2017), pp. 615–623.NDUSTRY-RELEVANT IMPLICIT LARGE-EDDY SIMULATION [53] R. C. Moura, J. Peir´o, and S. J. Sherwin , Implicit LES approaches via discontinuousGalerkin methods at very large Reynolds , in Direct and Large-Eddy Simulation XI, M. Sal-vetti, V. Armenio, J. Frohlich, B. Geurts, and H. Kuerten, eds., Springer, 2019, pp. 53–59.[54]
R. C. Moura, S. J. Sherwin, and J. Peir´o , Linear dispersion-diffusion analysis and its appli-cation to under-resolved turbulence simulations using discontinuous Galerkin spectral/ hp methods , Journal of Computational Physics, 298 (2015), pp. 695–710.[55] R. C. Moura, S. J. Sherwin, and J. Peir´o , Modified equation analysis for the discontinuousGalerkin formulation , in Spectral and High Order Methods for Partial Differential Equa-tions - ICOSAHOM 2014, R. M. Kirby, M. Berzins, and J. S. Hesthaven, eds., Springer,2015, pp. 375–383.[56]
R. C. Moura, S. J. Sherwin, and J. Peir´o , Eigensolution analysis of spectral/ hp continuousGalerkin approximations to advection-diffusion problems: insights into spectral vanishingviscosity , Journal of Computational Physics, 307 (2016), pp. 401–422.[57] D. Moxey, C. Cantwell, R. Kirby, and S. Sherwin , Optimising the performance of thespectral/ hp element method with collective linear algebra operations , Computer Methodsin Applied Mechanics and Engineering, 310 (2016), pp. 628–645.[58] D. Moxey, C. D. Cantwell, Y. Bao, A. Cassinelli, G. Castiglioni, S. Chun, E. Juda,E. Kazemi, K. Lackhove, J. Marcon, G. Mengaldo, D. Serson, M. Turner, H. Xu,J. Peir´o, R. M. Kirby, and S. J. Sherwin , Nektar++: enhancing the capability and ap-plication of high-fidelity spectral/hp element methods , Computer Physics Communications,249 (2020), p. 107110.[59]
D. Moxey, D. Ekelschot, ¨U.. Keskin, S. J. Sherwin, and J. Peir´o , A thermo-elastic analogyfor high-order curvilinear meshing with control of mesh validity and quality , Procediaengineering, 82 (2014), pp. 127–135.[60]
D. Moxey, M. D. Green, S. J. Sherwin, and J. Peir´o , An isoparametric approach to high-order curvilinear boundary-layer meshing , Computer Methods in Applied Mechanics andEngineering, 283 (2015), pp. 636–650.[61]
J. Nash , C1 isometric imbeddings , Annals of mathematics, (1954), pp. 383–396.[62]
Nektar++
OpenCascade
P.-O. Persson and J. Peraire , Curved mesh generation and mesh refinement using La-grangian solid mechanics , in 47th AIAA Aerospace Sciences Meeting including The NewHorizons Forum and Aerospace Exposition, 2009, p. 949.[65]
Pointwise
S. B. Pope , Ten questions concerning the large-eddy simulation of turbulent flows , New Journalof Physics, 6 (2004), p. 35.[67]
R. Poya, R. Sevilla, and A. J. Gil , A unified approach for a posteriori high-order curvedmesh generation using solid mechanics , Computational Mechanics, 58 (2016), pp. 457–490.[68]
M. Pratt , Introduction to ISO 10303 – the STEP standard for product data exchange , Journalof Computing and Information Science in Engineering, 1 (2001), pp. 102–103.[69]
O. Sahni, X. J. Luo, K. E. Jansen, and M. S. Shephard , Curved boundary layer meshingfor adaptive viscous flow simulations , Finite Elements in Analysis and Design, 46 (2010),pp. 132–139.[70]
D. Serson, J. R. Meneghini, and S. J. Sherwin , Direct numerical simulations of the flowaround wings with spanwise waviness , Journal of Fluid Mechanics, 826 (2017), pp. 714–731.[71]
M. S. Shephard, J. E. Flaherty, K. E. Jansen, X. Li, X. Luo, N. Chevaugeon, J.-F.Remacle, M. W. Beall, and R. M. O’Bara , Adaptive mesh generation for curved do-mains , Applied Numerical Mathematics, 52 (2005), pp. 251–271.[72]
S. J. Sherwin and M. Casarin , Low-energy basis preconditioning for elliptic substructuredsolvers based on unstructured spectral/ hp element discretization , Journal of ComputationalPhysics, 171 (2001), pp. 394–417.[73] S. J. Sherwin and J. Peir´o , Mesh generation in curvilinear domains using high-order ele-ments , International Journal for Numerical Methods in Engineering, 53 (2002), pp. 207–223.[74]
J. Slotnick, A. Khodadoust, J. Alonso, D. Darmofal, W. Gropp, E. Lurie, andD. Mavriplis , CFD vision 2030 study: a path to revolutionary computational aerosciences ,Technical Report CR 2014-218178, NASA, 2014.[75]
E. Tadmor , Convergence of spectral methods for nonlinear conservation laws , SIAM Journalon Numerical Analysis, 26 (1989), pp. 30–44.[76]
M. Taylor . personal communication, 2017.[77]
T. Toulorge, C. Geuzaine, J.-F. Remacle, and J. Lambrechts , Robust untangling of curvi-linear meshes , Journal of Computational Physics, 254 (2013), pp. 8–26. G. MENGALDO, ET AL.[78]
H. M. Tufo and P. F. Fischer , Fast parallel direct solvers for coarse grid problems , Journalof Parallel and Distributed Computing, 61 (2001), pp. 151–177.[79]
M. Turner , High-order mesh generation for CFD solvers , PhD thesis, Imperial College Lon-don, 2017.[80]
M. Turner, J. Peir´o, and D. Moxey , Curvilinear mesh generation using a variational frame-work , Computer-Aided Design, 103 (2018), pp. 73–91.[81]
B. C. Vermeire, S. Nadarajah, and P. G. Tucker , Implicit large eddy simulation usingthe high-order correction procedure via reconstruction scheme , International Journal forNumerical Methods in Fluids, 82 (2016), pp. 231–260.[82]
P. E. J. Vos, C. Eskilsson, A. Bolis, S. Chun, R. M. Kirby, and S. J. Sherwin , A genericframework for time-stepping partial differential equations (PDEs): general linear methods,object-oriented implementation and application to fluid problems , International Journal ofComputational Fluid Dynamics, 25 (2011), pp. 107–125.[83]
A. R. Winters, R. C. Moura, G. Mengaldo, G. J. Gassner, S. Walch, J. Peir´o, andS. J. Sherwin , A comparative study on polynomial dealiasing and split form discontinuousGalerkin schemes for under-resolved turbulence computations , Journal of ComputationalPhysics, 372 (2018), pp. 1–21.[84]
Z. Q. Xie, R. Sevilla, O. Hassan, and K. Morgan , The generation of arbitrary order curvedmeshes for 3D finite element analysis , Computational Mechanics, 51 (2013), pp. 361–374.[85]
C. Zhang, C. P. Bounds, L. Foster, and M. Uddin , Turbulence modeling effects on the CFDpredictions of flow over a detailed full-scale sedan vehicle , Fluids, 4 (2019), p. 148.[86]