Solving the incompressible Navier-Stokes equations is computationally heavy, with the pressure Poisson equation being the most time-consuming step. Iterative linear solvers are typically utilized to solve this equation. Since most solvers are iterative and rely on an initial guess, an opportunity emerges to use machine learning to improve this initial guess, such that fewer iterations are needed, consequently saving time.
A novel graph neural network is designed that employs a custom convolution algorithm, message-passing scheme, and pooling algorithm to maximize its performance. First, a convolution algorithm is proposed that uses interpolation to make a discrete (3x3) CNN kernel continuous. Then, instead of directly computing the kernel weight from the function, an integral over specified bounds is applied to account for the geometrical inhomogeneous distribution of the source nodes. The integral is embedded as the weighted sum of a vector containing learnable parameters, computed through a dot product with the edge attribute vectors. Next to the convolution operation, a message-passing scheme is designed that is compatible with the data format of the finite volume method whilst performing well in terms of the distance information can travel over the mesh. To conclude the design of the model, a custom pooling algorithm is designed that is equivalent to average pooling in CNNs.
A normalization procedure is established that ensures consistency in the model's magnitude. Notably, the ground truth output pressure is normalized using its standard deviation, which is unknown. To estimate this normalization factor, a correction model is established that uses the same convolution algorithm but employs an architecture inspired by classification CNNs.
The model's performance is evaluated based on the reduction in number of iterations required to reach convergence. This is done for both the Preconditioned Conjugate Gradient (PCG) solver and the multigrid Geometric Agglomerated Algebraic Multigrid (GAMG) solver
Across the various tests conducted, the number of iterations needed to reach convergence is reduced by approximately 40%, with the PCG solver performing slightly better than the GAMG solver. However, the PCG solver yields less consistent results, performing very well at samples that closely align with the training data, leading to a reduction of up to 60%. However, its performance drops significantly when tested on data that does not closely resemble the training data, sometimes even increasing the number of iterations. The GAMG solver demonstrates consistent performance, with almost no difference between the training and evaluation data.
In terms of generalization, the model demonstrates promising results, achieving similar performance across datasets with varying levels of complexity. This is interesting as the root mean square error differs significantly across the datasets and individual samples. This suggests that there is no direct relationship between the reduction in the number of iterations and the accuracy of the prediction. Furthermore, the model performs well on unseen meshes, showing that it can handle the unstructured nature of the meshes used throughout this research. This demonstrates the model's ability to be trained on a diverse dataset, after which it can be applied to unseen cases.

Among error estimation methods, adjoint-based approaches are considered the most accurate but have the highest computational cost. For unsteady non-linear problems such as the Navier-Stokes equations, substantial storage requirements arise, as the full primal solution must be stored to solve the adjoint problem. As a result, careful management of storage resources is essential. The purpose of this research was twofold: to develop a surrogate model for the primal solution and to compute an accurate adjoint-based error estimate with the developed surrogate primal.
This study compared three methodologies to create a surrogate model of the primal solution while reducing the storage requirements of unsteady adjoint-based error estimation: a Convolutional AutoEncoder (CAE), an Echo State Network (ESN) and a combination of the first two, referred to as CAE-ESN. A benchmark Proper Orthogonal Decomposition (POD) served as a baseline for comparison with the deep learning techniques. Three numerical test cases were analyzed, where the finite element method was used for spatial discretization and implemented with the \texttt{FEniCS} computational framework. The first test case involved a manufactured solution to verify the implemented solver and methodologies. The remaining test cases used a turbulent channel flow dataset to force the unsteady viscous Burgers' equations in 1D (wall-normal component) and 2D (spanwise and wall-normal components).
For the manufactured solution, the ESN and POD outperformed the remaining approaches for the lowest and highest spatial resolutions, respectively. The success of the ESN was linked to its training being a linear regression problem. As established in previous studies, the smooth nature of the solution rendered the POD optimal. For the 1D case, the CAE was optimal, particularly for lower spatial resolutions. This method offered equivalent compression ratios to the POD while being more efficient in terms of computational cost and accuracy. In contrast, the ESN-based methods failed to accurately capture the error estimate, as they were not able to accurately compute the primal residual. However, these methods offered a higher compression than other approaches, along with a decrease in accuracy. Moreover, the error indicators produced by the ESN-based methods continued to effectively pinpoint the elements required for mesh adaptation. In the 2D case, only the POD and CAE were investigated. The ESN was excluded due to the high-dimensional nature of the test case. The CAE-ESN was not applied because of the limited time interval, which provided insufficient data for training. The CAE again proved optimal due to its efficiency and higher compression capabilities than the POD. While both methods provided accurate error estimates and indicator fields, the CAE outperformed the POD due to its superior compression. The CAE was also able to compute the adjoint solution and primal residual more accurately than the POD for most spatial resolutions. This research highlighted the potential for the CAE to outperform more conventional methods, such as POD.

Goal-oriented reduced-order models (GOROM) have the potential to represent complex dynamics in an interpretable way. In the GOROM procedure, modes are found which when used in a Galerkin projection of the governing equations, return the desired goal function with minimum error. Determining such models leads to a high-dimensional optimisation problem, which may have multiple extrema, making purely local optimisation inapplicable.
In this thesis, a global optimisation algorithm was derived for GOROM. The optimisation response surface was approximated with a surrogate model, built using the Kriging technique. In each iteration, the surrogate model was improved by including the objective function of the three points with the highest expected improvements. This iterative process was stopped when the expected improvement fell below a threshold.
First, the performance of the algorithm was confirmed for a simple non-GOROM optimisation. Then the algorithm was applied to GOROM optimisation for the forced Burgers equation, initially using a complex forcing function giving multiple minima and later to a forcing function derived from a DNS of a turbulent channel flow. The initial forcing required on the order of tens of degrees of freedom to be optimised. The corresponding response surface proved similar to that of the simple problems, showing these optimisations to be tractable.
For the DNS forcing, hundreds of coefficients needed to be optimised. Using the developed algorithm a better optimal was found than with the local algorithm, although substantially more computational resources would be required to confirm global optimality.

This thesis presents a data-driven closure model for the variational multiscale method. The model is trained and tested in the context of a turbulent channel flow with Reτ = 180. Focus is on predicting the closure terms of the momentum equations. The model is trained using norms which only take into consideration local flow accuracy. For this reason long term stability is not addressed explicitly and is not the focus of this work. A convolutional neural network is chosen for its ability to efficiently make use of local spatial and temporal data. Different architectures and inputs are considered for the model. The performance is analyzed in an a priori sense by plotting correlations between the predicted closure terms and the exact closure terms. Using integrated forms of the momentum equation as input produces the highest a priori correlations. Taking multiple time steps as input also proves to be important. The model is then integrated into the simulation of a turbulent channel flow. The closure terms produced by the model are more accurate than those of a modern algebraic model. The flow field produced by the model is thus closer to that of the DNS than the flow field of the algebraic model.

An approach to data-driven closure modelling in the framework of variational multiscale method for large eddy simulation is presented. A turbulent channel flow at the friction Reynolds number of 180 is used as a case study. Challenges in the modelling linked to the continuity closure terms force the scope of the study to be narrowed to only momentum Navier-Stokes equations. By using integrated forms of the resolved flow solution to train a multi-layer perceptron closure model, high a priori correlations with true closure terms can be achieved. Validation of the data-driven closure models a posteriori shows that the models are suitable for the computation of the velocity field, but fail to obtain a physical pressure solution. Accumulation of the model error leads to a decay of the kinetic energy of the mean flow, while the resolved turbulent fluctuations become excessively energized. Novel techniques such as data augmentation, mini-batch sum loss and filtering of the closure terms are presented and evaluated via large eddy simulations.

Maneuvering or Stability and Control (S&C) characteristics of an aircraft constitute the most challenging and expensive phase of its design process. The S&C design phase extends into the development process, sometimes leading to unexpected aerodynamic issues. Thus, identifying aerodynamic issues and their consequences, e.g. transonic buffet, in the early stages of development, is crucial. Recent advancements in Computational Fluid Dynamics (CFD) have expanded its utilization. However, despite improvements in CFD, the comprehensive representation of dynamic effects for all possible maneuvers is still unfeasible due to the high computational expense.

Reduced Order Models (ROMs) have been combined with CFD data to predict an aircraft's dynamics in all possible maneuvers. ROMs enable the efficient utilization of high-fidelity CFD data, providing valuable insights into flight dynamic effects. This thesis project took place at the Netherlands Aerospace Center (NLR). The NLR in cooperation with TUDelft, has developed a ROM method for predicting unsteady aerodynamic loads of air vehicles. The current ROM approach combines the Proper Orthogonal Decomposition (POD) of pressure distribution with a Long Short-Term Memory (LSTM) type Neural Network (NN). So far, the POD-LSTM ROM method predicts the pressure distribution well in the incompressible flow regime. However, the increased number of spatial POD modes required to accurately represent the shock discontinuities in a pressure distribution poses challenges to POD-LSTM ROM. This leads to high computational costs, rendering the application of the POD-LSTM ROM in transonic flows impractical. Therefore, this thesis aims to set the foundation for expanding the POD-LSTM ROM for predicting the pressure distribution over sections of the DLR-F22 model in transonic conditions.

This research introduced a novel approach to address the increased number of spatial POD modes needed to approximate shock discontinuities in transonic flows. The enriched Proper Orthogonal Decomposition (ePOD) method introduces an enrichment basis into the standard truncated POD basis. The enrichment basis explicitly accounts for pressure discontinuities caused by shock waves, allowing the standard basis to focus on representing the remaining pressure distribution. The results confirm that the ePOD reduces the DoF required to approximate pressure distribution in transonic flows.

An LSTM neural network was utilized to forecast the time-dependent coefficients and parameters of the enriched reduced-order basis in unseen flow conditions. The results also showed that the ePOD reduced the complexity of the time-variant parameters of the reduced-order basis compared to the standard POD with the same number of degrees of freedom (DoF), facilitating more efficient training of the neural network.

Wings with leading-edge (LE) tubercles have gained increasing attention over the past decade. Despite their impressive aerodynamic performance, the underlying flow control mechanisms of tubercles remain controversial. In this thesis, both experimental and theoretical approaches are employed to investigate the flow patterns of a tubercled wing at pre-stall and post-stall angles of attack (AoAs).

In the experimental study, 2D Particle Image Velocimetry (PIV) was used to measure flow patterns at cross-flow planes along the chord. At a pre-stall AoA, high-vorticity regions generated by the tubercles appear in an alternating pattern near the LE. A quantitative comparison was conducted to examine the similarities between a tubercle and a delta wing. The results show that tubercles cannot be regarded as small delta wings in terms of vortex generation. The leading-edge vortex (LEV) sheets are convected downstream, where they interact with laminar separation bubbles (LSBs), creating complex flow patterns in the downstream regions. At a post-stall AoA, stall cells (SCs) appear along the span, with their formation dependent on both Reynolds number (Re) and tubercle amplitude. However, the spacing of SCs is relatively independent of AoA, Re, and amplitude, consistently ranging between 5 to 7 tubercle wavelengths.

In the theoretical study, the lifting line theory (LLT) approach was first used to predict the LEV strength but proved ineffective due to the absence of thickness effects. A subsequent analysis using the panel method in xflr5 showed that the Kutta condition should also be applied to the leading edge (LE) rather than only to the trailing edge (TE). Crow’s model was adapted by taking LEVs into consideration. However, a global description of the instability was not obtained due to difficulties in representing LEVs and related mathematical challenges.

This thesis contributes to a further understanding of the tubercle’s role in flow control. The LEVs generated by the tubercles are identified as key factors influencing flow evolution, yet these effects are not captured by LLT-based models or a conventional panel method. Future reduced-order models (ROMs) should account for the influence of LEVs to provide accurate representations of tubercled wing flow dynamics.

Unsteady numerical simulation has been proven to be an essential tool for research. The quality of the results can be improved by using mesh adaptation. Mesh adaptation uses error indicators to refine the mesh in regions with high errors. The error indicators used are output errors with the most accurate output error estimation method being adjoint-based error estimation. However, for this method, the primal solution needs to be stored, which is storage intensive, especially for large unsteady simulations. The method proposed in this thesis uses a neural network autoencoder to compress and reconstruct the primal solution. This solution is compared to a reconstructed solution using Proper Orthogonal Decomposition (POD).

The one-dimensional unsteady Burgers equation is used as validation for the methods using a manufactured solution while the lid-driven cavity flow is investigated using the proposed method. The manufactured solution of the one-dimensional Burgers case could be exactly reconstructed using two POD modes. For the autoencoder a small latent space was used. For low resolutions, the small latent space did not prove to be a problem as the primal and residual could be captured accurately. However, for higher resolutions, the reconstruction error of the autoencoder became dominant for the residuals and resulted in erroneous adjoint-based error estimates while the primal remained qualitatively similar.

For the lid-driven cavity flow, the POD was still able to capture the solution using a low number of modes due to the smoothness of the solution. This resulted in an unfair comparison between the POD and autoencoder reconstructed solutions. The reconstructed autoencoder error estimates for lower resolutions were more accurate due to the latent space being large enough to capture the residual of the discrete primal accurately enough. When moving to higher resolutions, the autoencoder was not able to reconstruct the residual accurately enough leading to erroneous error estimates. Therefore, the latent space of autoencoders should be sufficiently large in order to gain an accurate reconstruction of the residual. If the latent space is large enough, the error estimate is accurate and the local error estimates can be used as a first iteration error indicator for mesh refinement.

This thesis presents an energy-conservative data-driven approach in modelling the closure terms of the Navier-Stokes equations casted through the Variational Multiscale (VMS) framework. For context, the VMS framework is applied in designing stabilised finite element methods for multiscale phenomena in which stability is not guaranteed. Under standard applications of this framework with appropriate cubature rules and iteration routines, the only error source is model related in the assumptions for the fine scales, thereby influencing the closure terms of the coarse-scale equations. This issue can be rectified by predicting the perfect closure terms with the application of neural networks. To maintain energy conservation (as a generality, not as a guarantee), a constrained objective function was formulated, which was then translated into an unconstrained problem suitable for present day neural network training tools via the Augmented Lagrangian method. The features of the eight neural networks used were extracted from two spatiotemporal stencils. A priori evaluation of the fine-scale model demonstrated a tendency to overpredict closure terms for the Re_tau = 180 turbulent channel flow test case, yielding desired energy dissipative behaviour. This was true for both regimes of the flow in which the closure terms yield energy gain and energy loss via backscatter and cascade respectively. Unverified a posteriori simulation of the test case with the model displayed excessively high energy dissipation, in addition to high frequency oscillations, which were undampened by the time marching method used.

This research investigates the formation, growth, and size distribution of H2SO4 aerosols in aircraft engine plumes. It aims to enhance calculations by incorporating spatial variation and turbulence modeling. The study develops a toolchain using AER 3-D software, sets up a CFD model of an engine wake, and couples it with microphysics and advection schemes. The toolchain compares the AER 3-D microphysics model with state-of-the-art box models, showing its effectiveness in handling high H2SO4 concentrations and producing similar size distribution results. Using ANSYS, an axisymmetric engine wake field is calculated, revealing that condensation and nucleation occur in the injection area and wake boundary layer, while coagulation primarily occurs inside the wake. Turbulence influences aerosol diffusion and growth, resulting in size variations with distance from the wake center. The research emphasizes the importance of incorporating spatial variation and turbulence modeling for accurate aerosol predictions.

Stratospheric Aerosol Injection (SAI) is a geoengineering method to mitigate the effects of increased greenhouse gas concentrations in the Earth’s atmosphere, and to prevent further global warming. SAI does not reverse climate change, it merely counteracts its symptoms by offsetting the radiative forcing of greenhouse gases, and it should be accompanied by aggressive programmes of Carbon Dioxide Removal to cool down the climate. Operational studies suggest that the injection of condensable vapor from specialized high-altitude aircraft would be a reasonable option to form the desired aerosol particles. Yet, remarkably little is known about the growth evolution of aerosol inside a jet engine wake, especially in light of SAI where high initial concentrations of condensable vapor are injected into the exhaust stream. The lack of resolution in the flow field obscures a lot of the intricacies of the aerosol formation process, and raises serious doubts on the steerability of SAI, which is the capacity to spatially and temporally control its effects both in nature and scale. Without reassurance that the steerability requirement can be satisfied, the potential risks to the global ecosystem would be too high. In this respect, the study investigates the mechanisms within the near-field of a jet engine wake that lead to the creation of aerosols and subsequent growth, when accounting for local variations in temperature and relative humidity. This is done through a decoupled plume dispersion model that includes a flow solver which resolves the average velocity and turbulent intensity of the flow, a sectional chemistry module which includes all the relevant microphysical processes that affect aerosol growth, and a diffusion-advection model which calculates the displacement of aerosol, vapor, and chemiions in the flow field. In addition, several modeling errors are identified in classical thermodynamic approaches to high-density aerosol formation. Results demonstrate that the onset of particle formation in the plume is complex, and is heavily dependent on the injection concentration of sulfuric acid vapor and the mixing rate of the plume. Different aerosol particle sizes are formed depending on the location in the plume, leading to a nonuniform distribution of the volume-mean radius across the plume’s cross-sectional area. There is also a
strong indication that core particles experience preferential growth, which implies that the final particle distribution at the end of the plume’s lifetime might not be as uni-modal as self-limited theory predicts. The study concludes with a list of recommendations to solve the challenges that remain in order to truly
understand the early growth evolution of sulfate aerosol in a jet engine wake.

Hybrid RANS-LES methods have become a popular numerical approach for a wide variety of flows. This is due to dissatisfaction with the RANS modelling paradigm in separated flows along with the prohibitive computational cost of pure LES, especially in wall-bounded flows at high Reynolds numbers. However, these methods are susceptible to the grey area problem, where the modeling approach is neither RANS nor LES: rather it is a region with an ambiguous modeling approach. In zonal approaches that function as embedded wall-modeled LES (WMLES), the transition from RANS to LES can be accelerated by improving the synthetic turbulence and its injection into the flow. In this work, a systemic assessment of the two aspects of zonal grey area mitigation methods was carried out. The synthetic turbulence was generated by the synthetic turbulence generator (STG) and injected into the flow using two different forcing terms. To ensure accurate second-order statistics of synthetic turbulence, a priori estimations of the bias error associated with a specific realization of a random number set were implemented and used. This resulted in smaller deviations between the statistics of the synthetic turbulence and the target Reynolds stresses. Furthermore, a modified synthetic turbulence forcing that ensures more accurate estimation of the total shear stress in close proximity to the RANS-LES interface was proposed. Moreover, a dynamic forcing that selectively enhances the production of underestimated Reynolds stresses was implemented and evaluated. These aspects resulted in a faster transition from RANS to LES in terms of both skin friction coefficients and Reynolds stresses. In addition, the WMLES capabilities of the subgrid length scale Δ˜𝜔 together with the subgrid-scale 𝜎-model were explored. This work revealed that this combination is troublesome when used as embedded WMLES with synthetic turbulence, especially in stable flows. This is due to excessively decreased levels of eddy viscosity in the near-wall RANS region.

The aerodynamic model of a combat aircraft is essential for its success and competitiveness compared to other combat aircraft. This thesis aims to research the most optimal machine learning model to create an aerodynamic model of a combat aircraft. The very large but still sparse, highly nonlinear dataset forms a challenge for using specific machine learning models. Tree-based models, artificial neural networks (ANNs), and Bayesian neural networks (BNNs) have been identified to be individually capable of modelling the aerodynamics of combat aircraft. For ANNs, additional research was performed into genetic algorithms, as a robust hyperparameter optimisation method. Transfer learning, which reduced the training time by around 14%. Finally, adding gradient information of the training data as an additional input, reduced the mean squared error (MSE) by 7%. Scalable BNNs, which used Bayes-by-backprop, were developed to handle the large dataset. The uncertainties coming from the BNN were overconfident in the results compared to the tolerances of the dataset. The uncertainties showed only a weak correlation with the MSE of a given prediction.
Finally, to leverage each of the model's advantages, a stacked model was created which improves the predictive performance on average by 27% compared to the best base model in terms of the MSE on the test dataset. With the stacked model implemented, each of the aerodynamic coefficients could be predicted with over 97% of the predictions of the test data within the tolerance. This is an average increase of 0.5% for each of the aerodynamic coefficients compared to the best base models. Due to the strict regulations related to this aerodynamic model, the machine learning model that is created cannot replace the aerodynamic model at this time. However, the implementation of this machine learning model allows engineers to design the aerodynamic model faster, and with greater precision.

There is a need to examine the potential of alternative methods such as Variational Multiscale Method (VMM) for the ability to improve the efficiency of combustion simulations. However the computational expense of a LES run of a 3D turbulent non-premixed combustion problem is high due to the large variety of length and time scales in turbulent non-premixed combustion. Secondly, VMM is a sophisticated model and the models used to simulate turbulent non-premixed combustion are complex. And finally, the results of a full 3D turbulent non-premixed combustion simulation are difficult to interpret due to their complexity. Therefore the objectives as follows: Firstly to define a simplified model problem which has the important aspects of turbulent non-premixed combustion for comparing numerical methods. Secondly, to propose and to investigate the potential of a framework based on VMM formulation for non-premixed combustion computations. The results of the simplified model problem in a VMM framework are compared to the results when using a standard Smagorinsky model. Results show that in general the VMM model with the dynamic subscale assumption outperforms the Smagorinsky model when looking at the mean and RMS value of the velocity, mixture fraction and progress variable.

Reduced-order models (ROMs) are commonly used to reduce the complexity of large-scale fluid dynamics problems. The main technique to obtain their basis is proper orthogonal decomposition (POD), where an energy-based basis is obtained. Such a basis is derived without considering the governing equations. Moreover, if the flow dynamics are not energy-dominated, the resulting POD-ROM can be biased towards that quantity and not represent the case properly.

In this project, a constrained goal-oriented reduced-order modeling (GOROM) technique that optimizes the modes of POD-ROMs is applied to wall-bounded flows. Its semi-continuous formulation allows great flexibility in treating non-linear problems with discrete reference datasets and complex goal functionals, and since the model constraint is considered in the optimization, the output model has a more complete physical interpretation.

Two cases are considered. First, a 1D Burgers equation forced with data from a turbulent channel flow is studied. Different ROM sizes and goal functionals are tested and compared to the reference data statistics. After that, a 3D transitional channel is studied, with a reference dataset that is generated using an initial condition that experiences transient growth. In this case, the ROM size is fixed and the goal functionals used are related to this phenomenon.

The results for the Burgers case show the superiority of the GOROM compared to a POD-ROM. They also depend on the initial guess for the optimization modes, and it is shown how the choice of the goal functional influences the behavior of the solution. After that, the transitional channel highlights the influence of truncation on transient growth problems and shows the model constraint improves the growth prediction significantly. Finally, the use of dedicated goal functionals is shown to improve the prediction of relevant quantities such as the production of kinetic energy.

Vortex generators (VGs) are typically about two orders of magnitude smaller than their host component (such as an airplane wing). For this reason, conducting a fully-resolved RANS simulation to isolate their impact on the flow field is computationally expensive. This work presents a goal-oriented adjoint-based approach to model vortex generators using smooth source term distributions to make their simulations tractable.

Several approaches to modelling VGs have been proposed in the past. Phenomenological models such as BAY and jBAY form source term descriptions, while others directly model the vortex or its effects. Recently, a goal-oriented approach to model VGs was introduced. Florentie et al. optimized a set of discrete source terms in a cuboidal domain at the location of the VG, using a fully-resolved solution as a reference. The superior accuracy of the model motivated Dierickx et al. to develop a surrogate model of the optimized source term in order to extend its applicability to any mesh and account for the variability in boundary conditions (Inflow angle and Reynolds number). However, the reduced-order source term produced discretization errors due to the presence of steep gradients.

To overcome the shortcomings associated with the goal-oriented approach and its reduced-order modelling, in this work, a smooth-basis representation of a co-rotating VG array was constructed using a trust-region-based optimization algorithm. The descent direction was obtained from an adjoint system of equations formulated based on an unconstrained optimization problem involving an objective/goal functional and the RANS equations.

The developed framework produced source term distributions that were smooth enough to be resolved on a coarse mesh and accurate enough to generate a good representation of the VG’s effect on the flow field. Specifically, flow parameters such as the Reynolds stress, shape factor and boundary layer profiles corresponding to the fully-resolved simulation were accurately recreated on the coarse mesh. To account for the evolving inflow angle at the VG-surface junction during the flow field calculation, a surrogate model was constructed to vary the source term accordingly.

The developed smooth source term framework brings a new capability to aerodynamic design as, unlike the previous approaches, it can also be used to accurately model arbitrary flow control devices on relatively coarse meshes.

Reconstructed Discontinuous Galerkin (rDG) methods aim to provide a unified framework between Discontinuous Galerkin (DG) and finite volume (FV) methods. This unification leads to a new family of spatial discretization schemes from order three upwards. The first of these new schemes is the rDG(P1P2) method, which represents the solution on each element as linear functions while reconstructing quadratic contributions to compute the fluxes inside the element and over the faces. For the rDG(P1P2) method, two different reconstruction methods were implemented. The first of these reconstruction methods is a least-squares based reconstruction. For this reconstruction, an inverse distance weighting was introduced to improve the discretization error in anisotropic mesh regions, as well as an extended reconstruction stencil variant, which aims to stabilise the reconstruction on simplicial meshes. The inclusion of an inverse distance weighting was found to be beneficial for high Reynolds number flows on the example of the two-dimensional zero pressure gradient flat plate. As a second method, a variational reconstruction method was implemented. For the variational reconstruction rDG methods it was shown that they can offer significantly reduced discretization errors compared to DG methods for smooth flows. It was shown on the example of a method of manufactured solutions, that all implemented methods reach their designed order of accuracy and can provide lower spatial discretization errors than a DG method of a comparable order on regular and randomly perturbed hexahedral meshes as well as on tetrahedral meshes. The rDG methods was applied to several two and three-dimensional RANS test cases. For these test cases, a stronger influence of the Reynolds number on the discretization error of rDG methods was found compared to the weaker influence observed for DG methods. For all test-cases, it was shown that rDG methods converge faster on the same mesh, however, yield a higher absolute error, due to the lower number of degrees of freedom compared to native DG methods.

Structure-preserving or mimetic discretisations are a class of advanced discretisation techniques derived by employing concepts from differential geometry. Such techniques can attain specific conservation properties at the discrete level such as conservation of mass, kinetic energy, etc when applied to conservation laws. However, like traditional techniques, they are not entirely robust in specific multiscale cases such as under-resolved simulations which are often encountered in industrial applications. This ties into the concept of Large Eddy Simulations (LES) of which an enticing approach is the Variational Multiscale (VMS) method. Both Mimetic methods and VMS approaches have been extensively studied independently, however, their combination offers the potential of achieving a favorable robust solver capable of handling complicated industrial problems. Moreover, the extension of the Hybridised variant of the Mimetic methods towards advection-dominated problems is an interesting avenue yet to be explored. Therefore, this thesis focuses on the extension of the Hybridised Mimetic method and its combination with the VMS theory. The said combination is tested on simple linear equations, namely the advection-diffusion equation in both steady and unsteady cases. Additionally, the Hybridised Mimetic method is studied in more complex test cases such as the Burgers' equation and the incompressible Euler/Navier-Stokes equations.