Adaptive CFD schemes for aerospace propulsion

The flow fields which can be observed inside several components of aerospace propulsion systems are characterised by the presence of very localised phenomena (boundary layers, shock waves,...) which can deeply influence the performances of the system. In order to accurately evaluate these effects by means of Computational Fluid Dynamics (CFD) simulations, it is necessary to locally refine the computational mesh. In this way the degrees of freedom related to the discretisation are focused in the most interesting regions and the computational cost of the simulation remains acceptable. In the present work, a discontinuous Galerkin (DG) discretisation is used to numerically solve the equations which describe the flow field. The local nature of the DG reconstruction makes it possible to efficiently exploit several adaptive schemes in which the size of the elements (h-adaptivity) and the order of reconstruction (p-adaptivity) are locally changed. After a review of the main adaptation criteria, some examples related to compressible flows in turbomachinery are presented. An hybrid hp-adaptive algorithm is also proposed and compared with a standard h-adaptive scheme in terms of computational efficiency.


Introduction
Several components of aerospace propulsion systems are characterised by flow fields in which very localised phenomena can be observed. Shock waves, boundary layers, shear layers, turbulent coherent structures are some examples. These phenomena show strong gradients in space and they usually have a spatial extension which is reduced with respect to the size of the propulsion system. However, the performances of the propulsion system are strongly influenced by the effects of these structures. In order to accurately evaluate the flow field by means of Computational Fluid Dynamics (CFD) tools, it is necessary to refine sufficiently the computational mesh in the regions of interest. An efficient approach is based on automatic adaptive algorithms which can automatically refine the mesh only in the regions which require high resolution. This approach improves significantly the reliability of the simulation and there is a growing interest in the use of these techniques. Indeed, the design procedures adopted nowadays by the industries are often based on automatic optimisation algorithms which perform a systematic exploration of the design parameters. In this framework, it is fundamental to have tools which can automatically adapt the mesh to the several system configurations which are tested during the optimisation process. An overview of some adaptive techniques can be found in the results of the European Project ADIGMA (Adaptive Higher-Order Variational Methods for Aerodynamic Applications in Industry) [ Furthermore, the rapid development of High Performance Computing systems has recently introduced the possibility to perform large scale simulations based on massive parallelisation. While the classical industrial approach to fluid dynamics is based on Reynolds Averaged Navier-Stokes (RANS) simulations, Large Eddy Simulations (LES) are nowadays affordable [2]. The use of LES reduces the need of empirical turbulence models but requires a huge number of degrees of freedom. Furthermore, the fact that turbulent flow fields are unsteady makes it difficult to generate a good mesh a priori. For these reasons, several research efforts are nowadays devoted to the development of adaptive schemes.
In the present work, a Discontinuous Galerkin (DG) scheme is used to numerically solve the compressible Navier-Stokes equations. The main idea of the DG scheme is to describe the numerical solution inside each element as a linear combination of basis functions which do not have any continuity constraint at the interfaces between different elements. This feature makes the scheme extremely compact even when high-order reconstructions are employed. Furthermore, the local nature of the DG reconstruction is particularly useful when working with adaptive schemes. For example, non-conforming meshes (also known as hanging-nodes meshes) can be produced by adaptive algorithms based on element splitting and they can be easily managed by a DG scheme. This is a classical example of h-adaptivity, in which the size of the elements is locally changed in order to adapt the grid resolution to the local requirements. However, the DG discretisation offers also an alternative path to the development of adaptive schemes. Since the reconstruction order (p) inside each element is independent from the reconstruction order in the neighbouring elements, it is possible to locally change the order and to adapt the resolution level without changing the mesh. This approach, known as p-adaptivity, allows to arbitrarly choose the number of degrees of freedom (and so the order p) inside each element. In this paper an example of hybrid hp-adaptive algorithm is proposed and compared with a classical h-adaptive algorithm.

Space and time discretisation
The results presented in this work are obtained by a research code in which the governing equations are numerically solved by means of a semidiscrete approach. In particular, the equations are discretised in space by means of a discontinuous Galerkin method. The resulting ordinary differential equations are integrated in time by means of a fully implicit backward Euler scheme. This choice is related to the fact that only steady examples are presented in this work. The resulting linear system is solved by means of the restarted GMRES method with the ILU0 preconditioner provided by the free version of the PARALUTION library [3].

DG space discretisation
Consider a partial differential equation which describes a general conservation law: where u is the conservative variable, q is a source term and F is the flux vector (which can include both convective and diffusive terms). The solution u h can be described by a linear combination of N e basis functions: The Galerkin approach consists in setting to zero the projection of the residuals of the discretised equation on the space spanned by the basis functions.
The volume integrals which appear in Eq. 3 are computed on the element Ω e and the basis functions are defined only inside the element. For these reasons, there are not constraints related to the continuity of the solution between adjacent elements. The basis used in this work is a hierarchical modal basis obtained by applying the modified Gram-Schmidt orthonormalisation to a set of monomials defined in the physical space, according to Bassi et al. [4]. Local conservativity is guaranteed by using conservative numerical fluxes (F) across the interfaces. An upwind numerical flux based on the approach of [5] and [6] is used for the approximation of the convective fluxes. The diffusive fluxes are computed by means of a recovery based method [7]. The meshes are generated by the Gmsh tool [8] and then adapted in the CFD code.

Sensors for adaptive schemes
In order to adapt the discretisation according to the evolution of the solution it is necessary to use a sensor which can identifies the regions which need more resolution. In the framework of DG methods there are three main classes of sensors for adaptive schemes. The first class is based on the definition of an adjoint equation which relates the error on a target function with the local properties of the discretisation. In this way, it is possible to perform a goal-oriented optimisation of the mesh in order to minimise the discretisation error on a quantity of interest [9]. This approach can be extremely efficient in steady RANS simulations. The extension to unsteady problems is challenging because of large memory requirements and the coupling of space-time adaptation [10]. The second approach is based on feature-based indicators which can identify some characteristic features of the flow field by performing a check on some phisical variables. Even if featurebased indicators can appear as ad-hoc tools, they can be quite effective and extremely cheap to evaluate [11]. This makes them interesting for unsteady problems. Finally, it is possible to use the residuals of the discretisation to locally estimate the discretisation error on the conservative variables. This approach is not oriented to a particular target but aims at minimising all the sources of error in the field. This method is employed in this work. In particular the error indicator proposed by Leicht and Hartmann [12] is used: where h e is the characteristic element size, R and r are the volume and interface residuals respectively which are defined in [12].

An hybrid hp-adaptive scheme
The local nature of the DG reconstruction offers the possibility to easily implement both hadaptive and p-adaptive schemes. When a p-adaptive scheme is used, the reconstruction order is increased in the regions with the largest error sources: if the solution is sufficiently smooth, this approach leads to exponential convergence towards the exact solution. However, the regularity of the solution is often reduced by the presence of singularities in the flow field (for example shocks) or in the geometry (for example sharp corners). In these circumstances, high-order reconstructions are not particularly useful because they cannot get the theoretical accuracy order. For this reason, it could be better to switch to an h-adaptive algorithm which refines the mesh close to the non-smooth features. The convergence rate obtained by the hadaptive algorithm is not exponential but it is limited by the reconstruction order of the scheme and by the regularity of the solution. However, h-adaptivity is an effective way to keep under control the error related to a singularity in the flow field. The previous example suggests that an optimal apporach could be obtained by an hybrid algorithm which uses both h-adaptivity and p-adaptivity. In particular, the residual based error indicator described in the previous section is used to identify the elements which require more accuracy. Then it is necessary to decide if this element should be refined by splitting (h-adaptivity) or by increasing the reconstruction order (p-adaptivity). In order to make this choice it is possible to use an element wise smoothness sensor. Persson and Peraire [13] proposed a smoothness sensor for controlling artificial viscosity in shock capturing methods. Here, this sensor is applied to pressure in order to control the hybrid scheme: where P andP represent the pressure computed with the full DG expansion and with a truncated expansion in which the higher-order terms are neglected. The signal given by the sensor is related to the magnitude of the higher-order term in the expansion. These terms decay quickly on a smooth solution but they remain large in the presence of shock waves or other singularities. In the present work, the signal obtained by the smoothness indicator is compared with a threshold to decide whether the element should be refined by splitting or by p-refinement.

Transitional flow on the T106c cascade
In this Section the flow field in the T106c high-lift turbine cascade is studied. This test case shows several phenomena which are crucial for the design of low pressure gas turbines in turbofan engines. The working condition is described in [14] and is characterised by an inlet angle α i = 32.7 • and an isentropic exit Mach number M 2s = 0.65. The Reynolds number based on the axial chord and the isentropic exit conditions is Re 2s = 80000. The flow field is described by the RANS equations coupled with a three equations turbulence model [15]. A second order accurate DG scheme is used for these simulations. The chosen working condition is characterised by laminar separation on the suction side which is followed by transition to turbulence. In this test case it is fundamental to get enough resolution in the boundary layer and in the separated shear layer. The first requirement can be easily satisfied during the mesh generation process because it is sufficient to define the mesh element size as a function of the wall distance. However, the position of the shear layer is not know a priori and it is a result of the simulation. For this reason, it is useful to employ adaptive algorithms which can refine the mesh automatically during the simulation. In this test case an h-adaptive algorithm driven by the residual based error indicator of Section 3 is used. The original mesh and the adapted mesh after two refinement steps are shown in Fig. 1 and 2.
At each refinement step a fixed fraction of the elements with the largest error indicator are splitted. The entropy field (normalised with respect to the constant volume specific heat) and some streamlines obtained on the final adapted mesh are reported in Fig. 3.

Transonic flow on a NACA0012 airfoil
Shock waves can deeply influence the performances of compressors and high-pressure turbines in turbojet engines. In this test case, the transonic flow field around the NACA0012 airfoil is considered according to the working condition described in [16]. The fluid is assumed inviscid. The free stream Mach number and angle of attack are M ∞ = 0.8 and α ∞ = 1.25 • . A second order accurate (p=1) scheme is used to study the flow field and the mesh is adapted by using an h-adaptive algorithm based on the residual indicator described in Section 3. In Fig. 4    The same study is perfomed with a third order accurate DG scheme (p=2) and with the hybrid hp-adaptive algorithm described in Section 4. In Fig. 8 the error on the drag coefficient (computed with respect to a grid converged reference value) is reported as a function of the computational time for the different approaches. The plot shows that the implemented third order (p=2) scheme is not convenient with respect to the second order scheme (p=1) with the chosen adaptation strategy. This could be probably due to the presence of shock waves in the solution which limit the practical convergence order and reduce the benefits of the higher-order scheme. The hybrid scheme seems to be more efficient than the classical h-adaptive scheme for both reconstruction orders.

Conclusions
The main adaptation strategies for discontinuous Galerkin discretisations are discussed and their importance in the simulation of turbulent compressible flows in propulsion systems is highlighted.  Some applications of h-adaptive and hp-adaptive algorithms are presented. The proposed hybrid hp-adaptive algorithm is compared with an h-adaptive algorithm on a test case with shock waves. The hybrid algorithm performs better than the classical h-adaptive algorithm since it allows to get the same discretisation error with a lower computational cost. Future work will be devoted to dynamic load balancing to improve performances in adaptive parallel simulations.