High-order Computational Fluid Dynamics Research Articles

Evaluating performance portability of five shared-memory programming models using a high-order unstructured CFD solver

This paper presents implementing and optimizing a high-order unstructured computational fluid dynamics (CFD) solver using five shared-memory programming models: CUDA, OpenACC, OpenMP, Kokkos, and OP2. The study aims to evaluate the performance of these models on different hardware architectures, including NVIDIA GPUs, x86-based Intel/AMD, and Arm-based systems. The goal is to determine whether these models can provide developers with performance-portable solvers running efficiently on various architectures. The paper forms a more holistic view of a high-order solver across multiple platforms by visualizing performance portability (PP) and measuring productivity. It gives guidelines for translating existing codebases and their data structures to these models.

Verification of a High-Order FEM-based CFD Code using the Method of Manufactured Solutions

A high-order computational fluid dynamics (CFD) code capable of solving compressible turbulent flow problems was developed. The CFD code employs the Flowfield Dependent Variation (FDV) scheme implemented in a Finite Element Method (FEM) framework. The FDV scheme is basically derived from the Lax-Wendroff Scheme (LWS) involving the replacement of LWS’s explicit time derivatives with a weighted combination of explicit and implicit time derivatives. The code utilizes linear, quadratic and cubic isoparametric quadrilateral and hexahedral Lagrange finite elements with corresponding piecewise shape functions that have formal spatial accuracy of second-order, third-order and fourth-order, respectively. In this paper, the results of observed order-of-accuracy of the implemented FDV FEM-based CFD code involving grid and polynomial order refinements on uniform Cartesian grids are reported. The Method of Manufactured Solutions (MMS) is applied to governing 2-D Euler and Navier-Stokes equations for flow cases spanning both subsonic and supersonic flow regimes. Global discretization error analyses using discrete 𝐿2 norm show that the spatial order-of-accuracy of the FDV FEM-based CFD code converges to the shape function polynomial order plus one, in excellent agreement with theory. Uniquely, this procedure establishes the wider applicability of MMS in verifying the spatial accuracy of a high-order CFD code.

BROADCAST: A high-order compressible CFD toolbox for stability and sensitivity using Algorithmic Differentiation

The evolution of any complex dynamical system is described by its state derivative operators. However, the extraction of the exact N-order state derivative operators is often inaccurate and requires approximations. The open-source CFD (Computational Fluid Dynamics) code called BROADCAST discretises the compressible Navier-Stokes equations and then extracts the linearised N-derivative operators through Algorithmic Differentiation (AD) providing a toolbox for laminar flow dynamics. Furthermore, the gradients through adjoint derivation are extracted either by transposition of the linearised operator or through the backward mode of the AD tool. The software includes base-flow computation and linear global stability analysis via eigen-decomposition of the linearised operator or via singular value decomposition of the resolvent operator. Sensitivity tools as well as weakly nonlinear analysis complete the package. The numerical method for the spatial discretisation of the equations consists of a finite-difference high-order shock-capturing scheme applied within a finite volume framework on 2D curvilinear structured grids. The stability and sensitivity tools are demonstrated on two cases: a cylinder flow at low Mach number and a hypersonic boundary layer. Program summaryProgram Title: BROADCASTCPC Library link to program files:https://doi.org/10.17632/hcf893cr2n.1Developer's repository link:https://github.com/OneraHub/BroadcastLicensing provisions: Mozilla Public License 2.0Programming language: Python, FortranNature of problem: The extraction of the exact N-order state derivative operators useful to study complex dynamical systems is often inaccurate and requires approximations from the user.Solution method: This program discretises the compressible Navier-Stokes equations and then extracts the linearised state derivative operators through Algorithmic Differentiation providing a toolbox for laminar flow dynamics.

CFD Aided Ship Design and Helicopter Operation

In support of Canadian industrial and defence ship design and offshore helicopter operations, a series of Ship–Helicopter Operational Limits Analysis and Simulation (SHOLAS) projects are being conducted at the National Research Council Canada (NRC) in collaboration with Defence Research and Development Canada (DRDC). This study presents a brief overview of a Canadian in-house ship airwake simulation capability combining in-house high-fidelity wind-tunnel tests, full-scale sea trials, high-order computational fluid dynamics (CFD) tools, and realistic engineering-oriented flight simulators. This paper reports challenges and lessons learned during the course of the study, discusses the current capabilities and limitations of the CFD tools and the infrastructure required, and evaluates the gaps and barriers in industry adoption by focusing on how they could be overcome based on our current practice. After validating the CFD results of an updated version of a simplified frigate shape (SFS2) and the real-world Canadian Patrol Frigate (CPF), which are in reasonable agreement with the available in-house wind-tunnel and sea-trial data, the developed approach was recently applied to the design of an undisclosed Canadian ship. Among other applications, CFD airwake results were used with confidence as input to produce representative airwake features in industrial high-fidelity piloted flight simulators.

JAX-Fluids: A fully-differentiable high-order computational fluid dynamics solver for compressible two-phase flows

Physical systems are governed by partial differential equations (PDEs). The Navier-Stokes equations describe fluid flows and are representative of nonlinear physical systems with complex spatio-temporal interactions. Fluid flows are omnipresent in nature and engineering applications, and their accurate simulation is essential for providing insights into these processes. While PDEs are typically solved with numerical methods, the recent success of machine learning (ML) has shown that ML methods can provide novel avenues of finding solutions to PDEs. ML is becoming more and more present in computational fluid dynamics (CFD). However, up to this date, there does not exist a general-purpose ML-CFD package which provides 1) powerful state-of-the-art numerical methods, 2) seamless hybridization of ML with CFD, and 3) automatic differentiation (AD) capabilities. AD in particular is essential to ML-CFD research as it provides gradient information and enables optimization of preexisting and novel CFD models. In this work, we propose JAX-Fluids: a comprehensive fully-differentiable CFD Python solver for compressible two-phase flows. JAX-Fluids is intended for ML-supported CFD research. The framework allows the simulation of complex fluid dynamics with phenomena like three-dimensional turbulence, compressibility effects, and two-phase flows. Written entirely in JAX, it is straightforward to include existing ML models into the proposed framework. Furthermore, JAX-Fluids enables end-to-end optimization. I.e., ML models can be optimized with gradients that are backpropagated through the entire CFD algorithm, and therefore contain not only information of the underlying PDE but also of the applied numerical methods. We believe that a Python package like JAX-Fluids is crucial to facilitate research at the intersection of ML and CFD and may pave the way for an era of differentiable fluid dynamics. Program summaryProgram title: JAX-FluidsCPC Library link to program files:https://doi.org/10.17632/pzvkwn5s6p.1Developer's repository link:https://github.com/tumaer/JAXFLUIDSCode Ocean capsule:https://codeocean.com/capsule/6819679Licensing provisions: GNU GPLv3Programming language: PythonSupplementary material: Source code; Examples; Videos: Moving solid bodies, Taylor-Green vortex, Rising bubble, Shock-bubble interaction.Nature of problem: The compressible Navier-Stokes equations describe continuum-scale fluid flows. These flows often involve highly complex flow phenomena such as shocks, material interfaces, and turbulence. The intrinsic nonlinear dynamics render the numerical simulation of the these equations challenging. Machine learning provides novel avenues for describing partial differential equations. Machine learning models rely on gradient information provided by automatic differentiation and are often implemented in Python. In contrast, existing high-performance computational fluid dynamics codes are typically written in Fortran or C++ and do not offer inherent automatic differentiation capabilities. These discrepancies hinder the advance of machine-learning-supported computational fluid dynamics. Up to this day, a general-purpose fully-differentiable computational fluid dynamics solver for compressible two-phase flows is missing.Solution method: We introduce JAX-Fluids: a general-purpose three-dimensional fully-differentiable computational fluid dynamics solver for compressible two-phase flows. JAX-Fluids is a simulation framework intended for machine-learning-supported computational fluid dynamics research. Our framework is written entirely in JAX, a high-performance numerical computing library with automatic differentiation capabilities. We have used an object-oriented programming style and a modular design philosophy. This allows the straightforward exchange of numerical methods. We provide a wide variety of state-of-the-art high-order computational methods for compressible flows. The modularity of our framework additionally facilitates the integration of custom subroutines. We use the sharp-interface level-set method to model two-phase flows. The software package can easily be installed as a Python package. We have build the source code around the JAX NumPy API. This makes JAX-Fluids accessible and performant. JAX-Fluids runs on CPUs, GPUs, and TPUs. We use HDF5 in combination with XDMF for writing output quantities. The Python packages Haiku and Optax are used for implementation and training of machine learning methods.Additional comments including restrictions and unusual features: JAX-Fluids relies on open-source third-party Python libraries. These are automatically installed. In the current version, JAX-Fluids only runs on a single accelerator (CPU/GPU/TPU). Future versions will include support for parallel execution. JAX-Fluids has been tested on Linux and macOS operating systems.

High-Order CFD Solvers on Three-Dimensional Unstructured Meshes: Parallel Implementation of RKDG Method with WENO Limiter and Momentum Sources

In aerospace engineering, high-order computational fluid dynamics (CFD) solvers suitable for three-dimensional unstructured meshes are less developed than expected. The Runge–Kutta discontinuous Galerkin (RKDG) finite element method with compact weighted essentially non-oscillatory (WENO) limiters is one of the candidates, which might give high-order solutions on unstructured meshes. In this article, we provide an efficient parallel implementation of this method for simulating inviscid compressible flows. The implemented solvers are tested on lower-dimensional model problems and real three-dimensional engineering problems. Results of lower-dimensional problems validate the correctness and accuracy of these solvers. The capability of capturing complex flow structures even on coarse meshes is shown in the results of three-dimensional applications. For solving problems containing rotary wings, an unsteady momentum source model is incorporated into the solvers. Such a finite element/momentum source hybrid method eliminates the reliance on advanced mesh techniques, which might provide an efficient tool for studying rotorcraft aerodynamics.

Application of synthetic jet actuator to delay the flow seperation in symmetrical airfoil

The objective of this project is to delay the flow separation over the symmetrical airfoil by placing the synthetic jet actuator. A synthetic jet actuator is a device which moves a fluid in and outside a cavity through a small orifice and also it will maximize the aerodynamic performance of the airfoil. The aerodynamic characteristic of the airfoil have been compared in detail at different angles of attack, for the case with and without adoption of synthetic jet. This control device is having zero mass flux active flow control device. The external flow around the airfoil will be determined by using numerical approach. By changing the angle of attack we can analyse how the separation of flows are varied. Results of this study demonstrated the actuator effectiveness on overall aerodynamic performance and show consistent trends with high-order Computational Fluid Dynamics (CFD).

Thermal hydraulic modeling of solid fueled nuclear thermal propulsion reactors Part I: Development and verification

This paper focuses on the development and demonstration of a modern thermal–hydraulic engineering tool to design nuclear thermal propulsion reactors. Low enriched uranium nuclear thermal propulsion reactors are currently being invested to appease regulatory and proliferation concerns, and consequently rendering the previous highly enriched uranium designs inapplicable. Sparse experimental data from the previously designed highly enriched uranium designs may be non-ideal for benchmarking low-enriched uranium designs. Computational fluid dynamics codes provide an affordable and expeditious alternate route for the verification of thermal–hydraulic engineering tools. Verification of engineering tools is crucial as they will be implemented in coupled neutronic and thermal–hydraulic simulations to ascertain the thermal and safety margins of a given design. Extensive research was performed on the optimization and design of nuclear thermal propulsion reactors, but there is limited publicly available description of currently adopted thermal–hydraulic solution procedures used for the design of contemporary nuclear thermal propulsion reactors. More specifically, the modelling assumptions, geometric simplifications, correlations, and iterative procedures are not always well understood or documented. The purpose of this paper is to provide a clear explanation of our thermal hydraulic modelling approach, clarify the impacts of commonly used assumptions on engineering-based tools, and demonstrate the importance of thermal–hydraulic solutions on neutronics. These objectives are achieved through the newly developed ntpThermo computational tool. This tool is verified against the high-order computational fluid dynamics code OpenFoam, which serves as the reference solution following a successful validation against experimental data. The results indicate that the reduced order-modeling approach implemented in ntpThermo can provide valuable engineering design support and be used for coupled neutronic and thermal–hydraulic simulations.

A high-order partitioned fluid-structure interaction framework for vortex-induced vibration simulation

Vortex-induced vibration (VIV) is a highly nonlinear fluid-structure interaction (FSI) phenomenon, and its accurate prediction remains challenging. In this study, we developed a high-order partitioned FSI framework towards accurate simulation of the VIV phenomenon. The focus is on simulating the two-way coupled interaction between a flexibly mounted rigid body and its surrounding fluids. In the FSI framework, a high-order computational fluid dynamics method, i.e., flux reconstruction/correction procedure via reconstruction, is utilized to simulate the unsteady compressible Navier–Stokes equations in the arbitrary Lagrangian Eulerian format on unstructured moving grids. A linearized fluid-solid Riemann solver is employed to weakly couple fluid dynamics and rigid-body dynamics on the fluid-solid interface, which can effectively mitigate the so-called added-mass instability. Time marching is conducted through either the linearly implicit Rosenbrock–Wanner method or the explicit strong stability preserving Runge–Kutta method. The high-order FSI framework is verified with a challenging VIV phenomenon, i.e., VIV of a zero-mass cylinder in viscous flows, and has been used to study two VIV-related problems, namely, the nonlinear energy sink and oscillating foil based energy harvesting mechanism.

Gust mitigation with an oscillating airfoil at low Reynolds number

The encounter between micro-aerial vehicles (MAVs) and gusts is often detrimental and mitigating the effects of the gust is important for operating MAVs under severe environmental conditions. This study investigates the impact of vertical gusts on stationary and oscillating NACA0012 (National Advisory Committee for Aeronautics) airfoils at low Reynolds numbers using high-order computational fluid dynamics methods, and identifies key dynamics that dominate gust mitigation. The interaction of the gusts with the stationary airfoil generates large unsteady forces, which exceed the peak static lift coefficient. A simple pitch-down maneuver and oscillating airfoil motion were tested as methods for mitigating the effects of the gusts. A rapid and significant pitch-down maneuver is observed to inadvertently cause a stall event by exceeding the negative stall angle. A stepwise change in the angle of attack (AoA), as the gust develops, is shown to be more effective at mitigating the gust effect. However, this gust mitigation strategy is still not effective if the gust continues to grow in magnitude. Low amplitude wing oscillations were then tested as a novel method for gust mitigation. Increasing the reduced frequency of the oscillating airfoil is shown to dominate the gust and results in a predictable oscillatory lift and drag/thrust behavior. Results also show that this effect is relatively insensitive to variations in the Strouhal number. These results suggest there may be gust mitigation strategies leveraging oscillating wing behaviors on MAVs.

Python at Petascale With PyFR or: How I Learned to Stop Worrying and Love the Snake

In this work, we make the case for Python as a first-class language for high-performance computing (HPC) applications. Our rationale is based on three trends: an increasing desire for applications to be user friendly; the shift toward such applications becoming primarily consumers of third-party application programming interfaces; and growing use of code generation for performance sensitive kernels. In each instance, we argue why Python is better placed to ride these trends than traditional languages such as C or Fortran. Further, we will also address two additional concerns regarding the suitability of Python: long start-up times and limited support for multithreading. An as exemplar of Python in HPC, we consider the high-order computational fluid dynamics code PyFR. The novel architecture of PyFR is described including examples showcasing its domain-specific language.

High-order accurate direct numerical simulation of flow over a MTU-T161 low pressure turbine blade

Reynolds Averaged Navier-Stokes (RANS) simulations and wind tunnel testing have become the go-to tools for industrial design of Low-Pressure Turbine (LPT) blades. However, there is also an emerging interest in use of scale-resolving simulations, including Direct Numerical Simulations (DNS). These could generate insight and data to underpin development of improved RANS models for LPT design. Additionally, they could underpin a virtual LPT wind tunnel capability, that is cheaper, quicker, and more data-rich than experiments. The current study applies PyFR, a Python based Computational Fluid Dynamics (CFD) solver, to fifth-order accurate petascale DNS of compressible flow over a three-dimensional MTU-T161 LPT blade with diverging end walls at a Reynolds number of 200,000 on an unstructured mesh with over 11 billion degrees-of-freedom per equation. Various flow metrics, including isentropic Mach number distribution at mid-span, surface shear, and wake pressure losses are compared with available experimental data and found to be in agreement. Subsequently, a more detailed analysis of various flow features is presented. These include the separation/transition processes on both the suction and pressure sides of the blade, end-wall vortices, and wake evolution at various span-wise locations. The results, which constitute one of the largest and highest-fidelity CFD simulations ever conducted, demonstrate the potential of high-order accurate GPU-accelerated CFD as a tool for delivering industrial DNS of LPT blades.

Inline vector compression for computational physics

A novel inline data compression method is presented for single-precision vectors in three dimensions. The primary application of the method is for accelerating computational physics calculations where the throughput is bound by memory bandwidth. The scheme employs spherical polar coordinates, angle quantisation, and a bespoke floating-point representation of the magnitude to achieve a fixed compression ratio of 1.5. The anisotropy of this method is considered, along with companding and fractional splitting techniques to improve the efficiency of the representation. We evaluate the scheme numerically within the context of high-order computational fluid dynamics. For both the isentropic convecting vortex and the Taylor–Green vortex test cases, the results are found to be comparable to those without compression. Performance is evaluated for a vector addition kernel on an NVIDIA Titan V GPU; it is demonstrated that a speedup of 1.5 can be achieved.

Efficient Control-Oriented Coupled Electrochemical Thermal Modeling of Li-Ion Cells

Safe and effective exploitation of lithium ion (Li-ion) batteries requires advanced battery management systems. This article proposes a computationally efficient, control-oriented model of a Li-ion cell. The model describes the spatial nature of both the chemical species and temperature dynamics in a computationally efficient way. The method takes advantage of the algebraic structure that arises from the distributed nature of the model. We show that, by discretizing the model partial differential equations with a finite difference method, the coupling equations take a semiseparable structure for which an efficient algebra exists. This approach yields an efficient modeling tool that can be employed to design model-based estimation and control algorithms. The proposed model is validated against a high order computational fluid dynamics model showing accuracy and efficiency.

Enhanced Delayed Detached Eddy Simulation for Cavities and Labyrinth Seals

Abstract A range of popular hybrid Reynolds-averaged Navier–Stokes -large eddy simulation (RANS-LES) methods are tested for a cavity and two labyrinth seal geometries using an in-house high-order computational fluid dynamics (CFD) code and a commercial CFD code. The models include the standard Spalart–Allmaras (SA) and Menter shear stress transport (SST) versions of delayed detached eddy simulation (DDES) and the Menter scale adaptive simulation (SAS) model. A recently formulated, enhanced, variant of SA-DDES presented in the literature and a new variant using the Menter SST model are also investigated. The latter modify the original definition of the subgrid length scale used in standard DDES based on local vorticity and strain. For all geometries, the meshes are considered to be hybrid RANS-LES adequate. Very low levels of resolved turbulence and quasi-two-dimensional flow fields are observed for the standard DDES and SAS models even for the test cases here that contain obstacles, sharp edges, and swirling flow. Similar findings are observed for both the commercial and in-house high-order CFD codes. For the cavity simulations, when using standard DDES and SAS, there is a significant under prediction of turbulent statistics compared with experimental measurements. The enhanced versions of DDES, on the other hand, show a significant improvement and resolve turbulent content over a wide range of scales. Improved agreement with experimental measurements is also observed for profiles of the vertical velocity component. For the first labyrinth seal geometry swirl velocities are more accurately captured by the enhanced DDES versions. For the second labyrinth seal geometry, the mass flow coefficient prediction using the enhanced models is significantly improved (up to 30%). Standard, industrially available hybrid RANS-LES models, when applied to the present canonical cases can produce little to no resolved turbulent content. The standard SA- and Menter-based DDES models can yield lower levels of eddy viscosity when compared to equivilent steady RANS simulations which means that they are not operating as RANS or LES. It is recommended that hybrid RANS-LES models should be extensively tested for specific flow configurations and that special care is exercised by CFD practitioners when using many of the popular hybrid RANS-LES models that are currently available in commercial CFD packages.

Static aeroelastic stiffness optimization of a forward swept composite wing with CFD-corrected aero loads

This paper presents an aero load correction strategy applicable to the static aeroelastic optimization of composite wings. The optimization framework consists of a successive convex subproblem iteration procedure, in which a gradient-based optimizer consecutively solves a local approximation problem. Responses are approximated as a linear and/or reciprocal function of the laminate membrane and bending stiffness matrices. Together with the laminate thicknesses h, they constitute the design variables in the optimization process. Internally, the design space is transformed from stiffness matrices to lamination parameters, resulting in a continuous and convex optimization problem. Structural responses considered in the stiffness optimization are strength, local buckling and mass; aileron effectiveness, divergence, and twist constitute the aeroelastic responses. Steady aeroelastic loads are calculated with a doublet lattice method (DLM) embedded in the applied finite element solver, allowing for the generation of response sensitivities that incorporate the effects of displacement-dependent aeroelastic loads. To incorporate flow phenomena that cannot be reproduced with DLM, a higher order aerodynamic method is considered. The developed correction methods and their application are presented in this paper. The correction is twofold, first, aiming at a correction of DLM by means of camber and twist modifications applied directly to the doublet lattice mesh and second, by employing the capabilities of a higher order computational fluid dynamics (CFD) solver, like the DLR-based TAU code. To this end, DLM loads transferred to the structure are rectified by means of the supposedly superior CFD results. The aero load correction method is applied in the stiffness optimization of a forward swept wing. First, a trim application without structural optimization is discussed, to demonstrate the convergence behavior of the correction forces. The results of a wing skin mass minimization with balanced and unbalanced laminates are presented. In particular, the differences between optimizations with and without aero correction are highlighted. Eventually, a stacking sequence optimization based on the continuous optimization results is demonstrated.

A computational toolchain for the automatic generation of multiple Reduced-Order Models from CFD simulations

This paper describes the development of a systematic tool chain capable of automatically extracting accurate and efficient Reduced-Order Models (ROMs) from Computational Fluid Dynamics (CFD) simulations. These ROMs can then be used to support the design and operation of Near-Zero Energy Buildings (NZEB), with a higher accuracy than traditional zonal models but at a fraction of the computational cost of CFD. This study assesses the accuracy and time to solution of these ROMs when solved for appropriate Boundary Conditions (BCs), found in the built environment, in order to define the usability envelope of the automatically extracted ROMs. The parameters used in this study are inlet temperatures (K) and mass flow rates (kg/s). Results demonstrate that the absolute error can be maintained at under 0.5 K for changes in temperature of up to ±15 K, and under 0.25 K for changes in mass flow rates of up to ±45% of the original value. The results show that this method has the potential for applications in the built environment where the ROM accuracy and low computational cost can bridge a gap between low order RC models and high order CFD, further improving the energy efficiency in smart buildings.

High-order computational fluid dynamics simulations of a spinning golf ball

This paper presents the first high-order computational fluid dynamics (CFD) simulations of static and spinning golf balls at realistic flow conditions. The present results are shown to capture the complex fluid dynamics inside the dimples which lead to drag reduction versus a smooth sphere, and compare well to previous experimental and computational studies. The high-order flux reconstruction method has been paired with the artificial boundary overset method to enable simplified mesh generation and grid motion. The compressible Navier–Stokes equations are modeled using a scale-resolving large eddy simulation (LES) approach with no sub-grid models. The codes implementing these methods have been implemented for NVIDIA graphical processing units (GPUs), enabling large speedups over traditional computer hardware. The new method allows for the simulation of golf balls, and other objects at moderate Reynolds numbers, to be simulated in a matter of days on large computing clusters.

High-Order Flux Reconstruction on Stretched and Warped Meshes

High-order computational fluid dynamics is gathering a broadening interest as a future industrial tool, with one such approach being flux reconstruction (FR). However, due to the need to mesh complex geometries if FR is to displace current lower–order methods, FR will likely have to be applied to stretched and warped meshes. Therefore, it is proposed that the analytical and numerical behaviors of FR on deformed meshes for both the one-dimensional linear advection and the two-dimensional Euler equations are investigated. The analytical foundation of this work is based on a modified von Neumann analysis for linearly deformed grids, which is presented. The temporal stability limits for linear advection on such grids are also explored analytically and numerically, with Courant–Friedrichs–Lewy (CFL) limits set out for several Runge–Kutta schemes, with the primary trend being that contracting mesh regions give rise to higher CFL limits, whereas expansion leads to lower CFL limits. Lastly, the benchmarks of FR are compared to finite difference and finite volumes schemes, as are common in industry, with the comparison showing the increased wave propagating ability on warped and stretched meshes, and hence FR’s increased resilience to mesh deformation.

A family of high-order gas-kinetic schemes and its comparison with Riemann solver based high-order methods

Most high order computational fluid dynamics (CFD) methods for compressible flows are based on Riemann solver for the flux evaluation and Runge–Kutta (RK) time stepping technique for temporal accuracy. The advantage of this kind of space–time separation approach is the easy implementation and stability enhancement by introducing more middle stages. However, the nth-order time accuracy needs no less than n stages for the RK method, which can be very time and memory consuming due to the reconstruction at each stage for a high order method. On the other hand, the multi-stage multi-derivative (MSMD) method can be used to achieve the same order of time accuracy using less middle stages with the use of the time derivatives of the flux function. For traditional Riemann solver based CFD methods, the lack of time derivatives in the flux function prevents its direct implementation of the MSMD method. However, the gas kinetic scheme (GKS) provides such a time accurate evolution model. By combining the second-order or third-order GKS flux functions with the MSMD technique, a family of high order gas kinetic methods can be constructed. As an extension of the previous 2-stage 4th-order GKS, the 5th-order schemes with 2 and 3 stages will be developed in this paper. Based on the same 5th-order WENO reconstruction, the performance of gas kinetic schemes from the 2nd- to the 5th-order time accurate methods will be evaluated. The results show that the 5th-order scheme can achieve the theoretical order of accuracy for the Euler equations, and present accurate Navier–Stokes solutions as well due to the coupling of inviscid and viscous terms in the GKS formulation. In comparison with Riemann solver based 5th-order RK method, the high order GKS has advantages in terms of efficiency, accuracy, and robustness, for all test cases. The 4th- and 5th-order GKS have the same robustness as the 2nd-order scheme for the capturing of discontinuous solutions. The current high order MSMD GKS is a multi-dimensional scheme with incorporation of both normal and tangential spatial derivatives of flow variables at a cell interface in the flux evaluation. The scheme can be extended straightforwardly to viscous flow computation in unstructured mesh. It provides a promising direction for the development of high-order CFD methods for the computation of complex flows, such as turbulence and acoustics with shock interactions.