WP2

Main Achievement

WP2 - High accuracy numerical simulations and data management

Objectives: To improve existing high-order CFD numerical methods for the accurate prediction of unsteady flows on complex configurations. To improve their computational efficiency on HPC platforms by means of adaptation and implicit time-integration schemes. To develop new algorithms for data analysis and feature detection.

The short description of the activities carried out by the partners is given bellow. A more detailed description is provided in Deliverables D2.1 to D2.4.

T2.1 Development of high order schemes.

High-order numerical schemes are essential for solving differential equations accurately and efficiently. They provide greater precision with fewer grid points, reducing computational costs. These methods resolve small-scale features, such as turbulence and shock waves, with minimal numerical dissipation and dispersion errors, making them ideal for wave simulations. Despite requiring careful stability considerations, they offer better convergence and long-term integration accuracy. High-order schemes are widely used in computational fluid dynamics, electromagnetics, acoustics, and structural mechanics, improving simulations in aerospace, meteorology, and engineering.

In this task, various partners have collaborated to enhance the maturity of these schemes.

ESR 1 and ESR 16 investigated the applicability of mesh-free methods, particularly volume penalization or Immerse Boundary Methods (IBM). The Immersed Boundary Method (IBM) is a numerical technique used to simulate fluid flow around complex geometries without requiring a body-fitted mesh. Instead, the solid boundaries are represented within a fixed Cartesian grid, and their effects are imposed through additional forcing terms in the governing equations. This approach simplifies meshing, especially for moving or deformable bodies, and is widely used in computational fluid dynamics (CFD) for applications such as biofluid mechanics, aerodynamics, and multiphase flows. However, the limited performance of these approaches led to the numerical, physical, and mathematical development of a novel methodology that avoids high-order meshes while preserving accuracy.

ESR 1 and ESR 16 developed a IBM plugin into UPM HORSES3D high order solver. They compared the IBM results with body-fitted mesh simulations and reference data from the literature. This benchmarking phase was crucial in evaluating the accuracy and limitations of the IBM implementation, particularly regarding boundary layer effects and pressure distributions.

Building on these initial tests, they extended the study to a more complex three-dimensional aircraft wing in a landing configuration. A hybrid meshing strategy was adopted: the main wing was discretized using a conventional body-fitted mesh, while the trailing-edge flap was modeled with the IBM approach. The 3D meshing for this study was performed entirely in-house. The goal was to assess the feasibility and accuracy of IBM for geometrically complex components, particularly in scenarios where generating a high-quality body-fitted mesh is challenging.

The results from this hybrid approach were compared against a fully body-fitted mesh simulation to evaluate discrepancies in flow characteristics, aerodynamic forces, and computational efficiency.

As part of advancements in high-order methods, ESR4 has developed a proof-of-concept workflow for hrp-AMR (Adaptive Mesh Refinement) to address compressible and incompressible steady and unsteady flow simulations. ESR4 also explored various AMR techniques, including local solver expansion polynomial order (p-type), node movement (r-type), and element splitting (h-type).

A novel h-AMR isoparametric high-order splitting, swapping, and smoothing method was implemented in the open-source tool NekMesh, part of the Nektar++ meshing utility, enabling hrp-AMR. These new 2D and 3D mesh adaptation techniques allowed ESR4 to further develop h-r-p adaptation, particularly beneficial for compressible flows. In this methodology, hr-adaptation was used to capture discontinuous features such as shocks, while a variable p-order was applied in smooth regions.

The newly developed hrp-AMR approach was tested for improved efficiency and accuracy in capturing flows with complex shock physics. A key achievement was the hrp-AMR combination, which resulted in a 12-fold reduction in mesh count, outperforming conventional h-AMR techniques.

An example of this adaptation for NACA 0012 transonic flow at M = 0.8 is shown in Figure 2, while a comparison of hrp-adaptation versus standard h or hr-adaptation in terms of performance and efficiency can be seen in the table below

Figure 2 Example of hrp-AMR mesh with the shock sensor for a transonic case at M=0.8.

Additionally, ESR4 identified oscillations in high-order simulations caused by mesh geometrical accuracy, notably affecting shear stresses in TC1. A 2D study using TC1 Imperial Front Wing and NACA 0012 examined the required geometrical accuracy and polynomial order. Euler compressible flow at Mach 0.8, AoA = 1.25°, and Cp-distributions were analyzed for p1-p4 meshes. Results showed P1 meshes were inadequate, P2 meshes oscillated, and P4 meshes closely matched expected Cp-distributions. Findings motivated improvements in high-order mesh generation for WP2 – TC1 and NekMesh.

Figure 3 Instantaneous Cp-distributions on a transonic NACA0012 at M=0.8 and AoA=1.25° simulated with the same DG-Euler solver showing the oscillatory behavior of high order.

ESR7 developed and implemented a high-order transition model using the discontinuous Galerkin method, specifically designed for transonic airfoil flows. The model has been fully analytically linearized and integrated into the Aghora DG CFD solver framework. It has been validated across various flow problems, including zero-pressure-gradient flat plates and both steady and unsteady subsonic and transonic airfoil flows. The approach shows strong convergence properties and demonstrates favorable agreement with both experimental and computational results

ESR11 focused on demonstrating the feasibility of running the Fidelity High-Order Solver for automotive simulations, including the Imperial Front Wing (TC1) and the Windsor Body. The second objective was to establish meshing best practices for these test cases to showcase potential advantages over LES. However, due to internal management decisions at Cadence, the High-Order Solver project was discontinued, leading ESR11 to shift focus to a Finite Volume Wall-Modeled LES solver (CharLES)

ESR12 investigated different stabilization methods for high-order simulations during her secondment at UPM. She applied Gradient Jump Penalization and Spectral Vanishing Viscosity to simulate canonical cases, including 2D cylinder flow, 3D turbulent channel flow, and the Taylor-Green vortex, using Nektar++ with continuous approximations. ESR12 then compared her results with those from HORSES3D, which employs the same stabilization methods using a Discontinuous Galerkin projection.

ESR13 developed a novel preconditioning technique—the Lower-Order Refined (LOR) Preconditioner—to enhance the efficiency of solving the pressure Poisson system in high-order incompressible flow simulations. This approach uses a low-order (P=1) finite element discretization, which is spectrally equivalent to the high-order system. Key advantages include low-cost operator evaluations, constant memory per degree of freedom, robustness for complex geometries, and a controlled iterative condition number. The preconditioned system is solved using standard AMG routines, making it both easy to implement and highly efficient for large-scale industrial simulations.

Additionally, ESR13 played a crucial role in troubleshooting, performance testing, and code development as the first industrial user. He conducted scalability tests, implemented new functionalities, and improved post-processing capabilities, including large-scale I/O handling with HDF5 and parallel VTU outputs for HPC clusters.

T2.2 Temporal discretization in high order schemes.

Temporal discretization is vital in high-order schemes for solving time-dependent PDEs, ensuring accuracy, stability, and efficiency. Proper temporal discretization preserves solution quality in complex problems like shocks or multi-physics interactions, making it essential for robust, reliable simulations in computational fluid dynamics, and compatibility with high-order spatial discretization prevents accuracy bottlenecks.

ESR3 focused on developing a posteriori time error indicator for variable time step sizes within the BDF (Backward Differentiation Formula) family of methods. The proposed developments are also applicable to more complex time integration schemes. The main idea behind designing this a posteriori error indicator is that a better approximation of the time derivative can be achieved using the BDF scheme of order k+1 (BDF-(k+1)), even when applied to the sequence obtained with BDF-k. This approach provides an estimate of the error associated with the approximation of the time derivative.

Additionally, in the fractional step scheme ESR3 is using, a time error arises due to the splitting of the problem. ESR3 estimates this error by making it proportional to the difference between the intermediate velocity and the velocity at the end of the time step.

Figure 5 presents the convergence plots in which the fractional step scheme is employed for isentropic flows. The plots show the error predicted by our a posteriori error indicator and the true error. If the splitting error is not considered, the prediction (displayed in red) is poor, but when the estimate of the splitting error is introduced, the new prediction (displayed in dark blue) is very satisfactory.

ESR5 implemented an implicit temporal discretization for the incompressible Navier-Stokes equations. The advantage of an implicit scheme is that, while the computational cost per time step increases, the overall time-to-solution decreases due to increased numerical stability, which allows for larger time steps.

The time-stepping algorithm used is a fractional step scheme, known as velocity correction. This scheme decouples the pressure and velocity variables and employs a linearization of the advection term to ensure unconditional stability. We compared the performance of the implicit velocity correction scheme with that of an explicit scheme and evaluated the potential speed-up offered by the implicit time-stepping approach. The measurements were conducted using the Imperial Front Wing (IFW) geometry derived from TC1 (see deliverable D2.1 for more details).

ESR11 implemented a 1D code to solve simple problems like the advection equation, heat equation, and Euler equation, while also incorporating key components of Fidelity’s high-order solver, such as flux reconstruction discretization, dual time-stepping, and polynomial multigrid. The code supports various high-order temporal discretizations, including BDF and Adam-Moulton schemes, as well as explicit or implicit Runge-Kutta schemes. For many of these schemes, the theoretical accuracy was verified (see D2.3). However, for CharLES, we found that the temporal error was negligible compared to the spatial (meshing) error. For Wall-Resolved LES, which requires small time steps, using multigrid was found to be counterproductive. Multigrid is designed for better convergence of low-frequency modes, but with semi-implicit time stepping (CharLES), low frequencies remain almost unchanged between time steps, making multigrid unnecessary. Multigrid is more beneficial for fully implicit schemes.

ESR12 focused on the “Sub-Stepping” scheme, an implicit time-stepping method based on the semi-Lagrangian formulation in Nektar++. She first validated the scheme’s benefits on canonical flow problems and then applied it to industrial cases like the Extruded IFW and IFW. After identifying bottlenecks, the scheme’s computational performance was enhanced by using a matrix-free format. Parametric studies helped identify optimal configurations for industrial simulations. The Extruded IFW simulations were compared with those using ESR5’s implicit scheme, leading to a journal publication. Both schemes allowed for larger time steps, reducing computational time. While this resulted in reduced accuracy, the error in the integral forces was under 1% for time steps up to 20 times larger than the explicit time-integration limit, which is acceptable for industrial applications.

T2.3 Detailed simulation and control of TC1

The description here included is complementary with Task 1.2 in WP1.

ESR1 simulated a multielement configuration using both the standard volume penalization method (IBM) and a novel methodology that preserves high-order accuracy using the immersed boundary method. As observed in the figure, the latter provided the best results and was able to capture the underlying physics, particularly the laminar separation bubble

Figure 6 Comparison between the standard IBM (only first order), left, and the new implementation of high order IBM, right.

ESR3 initially focused on reproducing the simulation results of the front wing of TC1 in FEMUSS (a Finite Element flow solver, an in-house code developed at CIMNE). After confirming the results, the effort shifted towards testing the performance of the proposed a-posteriori time error indicator using an industrial-scale geometry with a fully turbulent flow regime.

For this test case, five different time step sizes, ranging from 0.0002 to 0.0016 (corresponding to a CFL number from 4.40 to 37.72), were selected and simulated over a time interval from 0 to 1 second, using the constant time step size BDF2 integration scheme. The simulations used the same fully developed initial condition.

Figure 7 Left: Close-up view of the velocity magnitude initial condition around the IFW. Right: illustrates the estimated time error throughout the simulation interval using the proposed time error indicator across five different simulations, each with varying constant time step sizes.

ESR4 applied h-adaptation in unsteady incompressible simulations, such as the extruded Imperial Front Wing (TC1), to improve the identification of flow transition to turbulence and the resolution of wake vortices. The study focused on feature-based error indicators combined with h-adaptation based on the mean Q-criterion. Starting with a coarse mesh, four adaptation cycles progressively refined the mesh around key regions, enhancing accuracy near the transition point and the suction side, while maintaining efficiency on the pressure side. The refined trailing edge mesh also captured flow sensitivity and propagating eddies. The same approach is planned for the full IFW TC1 simulation.

Figure 8 An example of the Q-criterion on a simulation of the extruded IFW TC1 and the final mesh after 4 adaptation cycles for capturing transition behavior of the eIFW TC1.

ESR5 collaborated with ESR13 on a detailed simulation of the front wing of TC1. The focus of this effort was to test time-stepping algorithms using an industrial-scale geometry with turbulent flow phenomena. The Imperial Front Wing (IFW) is an ideal test case, as it features challenging phenomena, such as laminar-to-turbulent transition in the boundary layer flow around the three wing profiles (see D2.1 for a detailed description of the geometry and flow physics discussion).

The computational mesh and simulation setup were validated against literature using data from a previous study. The comparison in Figure 2 shows the time-averaged skin friction coefficient (Cf) from the current effort and compares it with results from Slaughter et al. (2023). The comparison demonstrates good agreement for all three wing profiles. The transition to turbulence on the main plane wing around x/c = 0.5 and on the first flap around x/c = 1.25 is well captured. Additionally, the bypass transition process on the second flap, as described in Slaughter et al. (2023), was successfully identified.

ESR11 focused on establishing best practices for meshing the Extruded Imperial Front Wing using the CharLES solver, now used in automotive simulations. These practices were validated on the Imperial Front Wing (TC1) and compared to experimental data and CFD results from ESR12 and ESR13. During a secondment at McLaren, a rotating wheel was added. The impact of simulation timestep on flow accuracy was studied, revealing that larger timesteps reduced lift coefficient and altered vortex interactions. A method using implicit time-stepping for rapid flow transition was also explored.

T2.4 Algorithms for mode decomposition and data analysis.

Algorithms for mode decomposition and data analysis are crucial for extracting meaningful patterns and structures from complex datasets, especially in fluid dynamics, climate modeling, and signal processing. Techniques like Proper Orthogonal Decomposition (POD), Dynamic Mode Decomposition (DMD), and Singular Value Decomposition (SVD) allow the identification of dominant modes in the data, reducing dimensionality while preserving key features. These methods facilitate the analysis of dynamic systems by isolating temporal and spatial modes, aiding in the understanding of flow behaviors, oscillations, and turbulence. Mode decomposition plays a significant role in system prediction, control, and optimization in various engineering fields.

In this task:

ESR12 used statistical methods to assess convergence over time in simulations of industrial geometries, such as the Extruded IFW and IFW. Along with ESR13, they developed a Python library for post-processing integral forces over time, tracking the evolution of fields at probe points or 2D planes. This involved calculating time statistics, creating animations to visualize flow evolution, and transforming data into the frequency domain to identify dominant frequencies. They also utilized libraries with mode decomposition algorithms like SVD, POD, and DMD, specifically pyDMD and pyLOM.

ESR2 worked on two improvements to Variational Mode Decomposition (VMD). The first modification introduces objectives within the optimization problem to promote orthogonality between dominant modes. This is achieved by projecting components into other modes using coherence, combined with Wiener filtering, which improves convergence and uniqueness, particularly in over-segmentation conditions. This new method is called Orthogonalized Variational Mode Decomposition (OVMD).

The second improvement involves estimating an appropriate value for the weight parameter α, which influences the frequency sidebands of modes. This is done by identifying a representative bandwidth B from the input signal spectrum, and calculating αk as inversely proportional to the product of B and the mode’s central frequency. The new approach shows a stronger positive impact when combined with mode orthogonalization.

OVMD provides more accurate decomposition than VMD for broadband synthetic signals, especially under over-segmentation, and captures main modes more effectively, leaving extra modes with minimal amplitude or noise. Sensitivity assessments show that OVMD yields consistent results with various values of the inverse proportionality factor, noise levels, and extreme over-segmentation.

To obtain dominant features of the flow, ESR7 performed a global linear stability analysis for flow at Reynolds number Re = 5·10⁵ in the vicinity of transonic buffet. The eigenvalue spectrum was analyzed for increasing Mach numbers starting from M = 0.71. For M < 0.715, the flow is globally stable, but at M = 0.71, a low-frequency eigenvalue departs from the spectrum, indicating instability. The flow becomes globally unstable at M = 0.715 due to a Hopf bifurcation. The eigenvalue’s growth rate peaks at M = 0.75 and then decreases, with the flow stabilizing again at M = 0.8. The analysis aligns well with URANS and LES results, confirming that the unstable eigenvalue is responsible for the transonic buffet oscillations