Higher-order quadrature-based moment methods for kinetic equations

Kinetic equations containing terms for spatial transport, body forces, and particle–particle collisions occur in many applications (e.g., rarefied gases, dilute granular gases, fluid-particle flows). The direct numerical solution of the kinetic equation is usually intractable due to the large number of independent variables. A useful alternative is to reformulate the problem in terms of the moments of the velocity distribution function. Closure of the moment equations is challenging for flows sufficiently far away from the Maxwellian limit. In previous work, a quadrature-based third-order moment closure was derived for approximating solutions to the kinetic equation for arbitrary Knudsen number. A key component of quadrature-based closures is the moment-inversion algorithm used to find the non-negative weights and velocity abscissas. Here, a robust inversion procedure is proposed for three-component velocity moments up to ninth order. By reconstructing the velocity distribution function, the spatial fluxes in the moment equations are treated using a kinetic-based finite-volume solver. Because the quadrature-based moment method employs the moment transport equations directly instead of a discretized form of the kinetic equation, the mass, momentum and energy are conserved for arbitrary Knudsen and Mach numbers. The computational algorithm is tested for the Riemann shock problem and, for increasing Knudsen numbers (i.e. larger deviations from the Maxwellian limit), the accuracy of the moment closure is shown to be determined by the discrete representation of the spatial fluxes.     

Conditional quadrature method of moments for kinetic equations

Kinetic equations arise in a wide variety of physical systems and efficient numerical methods are needed for their solution. Moment methods are an important class of approximate models derived from kinetic equations, but require closure to truncate the moment set. In quadrature-based moment methods (QBMM), closure is achieved by inverting a finite set of moments to reconstruct a point distribution from which all unclosed moments (e.g. spatial fluxes) can be related to the finite moment set. In this work, a novel moment-inversion algorithm, based on 1-D adaptive quadrature of conditional velocity moments, is introduced and shown to always yield realizable distribution functions (i.e. non-negative quadrature weights). This conditional quadrature method of moments (CQMOM) can be used to compute exact N-point quadratures for multi-valued solutions (also known as the multi-variate truncated moment problem), and provides optimal approximations of continuous distributions. In order to control numerical errors arising in volume averaging and spatial transport, an adaptive 1-D quadrature algorithm is formulated for use with CQMOM. The use of adaptive CQMOM in the context of QBMM for the solution of kinetic equations is illustrated by applying it to problems involving particle trajectory crossing (i.e. collision-less systems), elastic and inelastic particle–particle collisions, and external forces (i.e. fluid drag).     

A high-order gas-kinetic Navier–Stokes flow solver

The foundation for the development of modern compressible flow solver is based on the Riemann solution of the inviscid Euler equations. The high-order schemes are basically related to high-order spatial interpolation or reconstruction. In order to overcome the low-order wave interaction mechanism due to the Riemann solution, the temporal accuracy of the scheme can be improved through the Runge–Kutta method, where the dynamic deficiencies in the first-order Riemann solution is alleviated through the sub-step spatial reconstruction in the Runge–Kutta process. The close coupling between the spatial and temporal evolution in the original nonlinear governing equations seems weakened due to its spatial and temporal decoupling. Many recently developed high-order methods require a Navier–Stokes flux function under piece-wise discontinuous high-order initial reconstruction. However, the piece-wise discontinuous initial data and the hyperbolic-parabolic nature of the Navier–Stokes equations seem inconsistent mathematically, such as the divergence of the viscous and heat conducting terms due to initial discontinuity. In this paper, based on the Boltzmann equation, we are going to present a time-dependent flux function from a high-order discontinuous reconstruction. The theoretical basis for such an approach is due to the fact that the Boltzmann equation has no specific requirement on the smoothness of the initial data and the kinetic equation has the mechanism to construct a dissipative wave structure starting from an initially discontinuous flow condition on a time scale being larger than the particle collision time. The current high-order flux evaluation method is an extension of the second-order gas-kinetic BGK scheme for the Navier–Stokes equations (BGK-NS). The novelty for the easy extension from a second-order to a higher order is due to the simple particle transport and collision mechanism on the microscopic level. This paper will present a hierarchy to construct such a high-order method. The necessity to couple spatial and temporal evolution nonlinearly in the flux evaluation can be clearly observed through the numerical performance of the scheme for the viscous flow computations.     

Optimal Runge–Kutta smoothers for the p-multigrid discontinuous Galerkin solution of the 1D Euler equations

This work presents a family of original Runge–Kutta methods specifically designed to be effective relaxation schemes in the numerical solution of the steady state solution of purely advective problems with a high-order accurate discontinuous Galerkin space discretization and a p-multigrid solution algorithm. The design criterion for the construction of the Runge–Kutta methods here developed is different form the one traditionally used to derive optimal Runge–Kutta smoothers for the h-multigrid algorithm, which are designed in order to provide a uniform damping of the error modes in the high-frequency range only. The method here proposed is instead designed in order to provide a variable amount of damping of the error modes over the entire frequency spectrum. The performance of the proposed schemes is assessed in the solution of the steady state quasi one-dimensional Euler equations for two test cases of increasing difficulty. Some preliminary results showing the performance for multidimensional applications are also presented.     

Nonuniform time-step Runge–Kutta discontinuous Galerkin method for Computational Aeroacoustics

With many superior features, Runge–Kutta discontinuous Galerkin method (RKDG), which adopts Discontinuous Galerkin method (DG) for space discretization and Runge–Kutta method (RK) for time integration, has been an attractive alternative to the finite difference based high-order Computational Aeroacoustics (CAA) approaches. However, when it comes to complex physical problems, especially the ones involving irregular geometries, the time step size of an explicit RK scheme is limited by the smallest grid size in the computational domain, demanding a high computational cost for obtaining time accurate numerical solutions in CAA. For computational efficiency, high-order RK method with nonuniform time step sizes on nonuniform meshes is developed in this paper. In order to ensure correct communication of solutions on the interfaces of grids with different time step sizes, the values at intermediate-stages of the Runge–Kutta time integration on the elements neighboring such interfaces are coupled with minimal dissipation and dispersion errors. Based upon the general form of an explicit p-stage RK scheme, a linear coupling procedure is proposed, with details on the coefficient matrices and execution steps at common time-levels and intermediate time-levels. Applications of the coupling procedures to Runge–Kutta schemes frequently used in simulation of fluid flow and acoustics are given, including the third-order TVD scheme, and low-storage low dissipation and low dispersion (LDDRK) schemes. In addition, an analysis on the stability of coupling procedures on a nonuniform grid is carried out. For validation, numerical experiments on one-dimensional and two-dimensional problems are presented to illustrate the stability and accuracy of proposed nonuniform time-step RKDG scheme, as well as the computational benefits it brings. Application to a one-dimensional nonlinear problem is also investigated.     

High-order finite-volume methods for the shallow-water equations on the sphere

This paper presents a third-order and fourth-order finite-volume method for solving the shallow-water equations on a non-orthogonal equiangular cubed-sphere grid. Such a grid is built upon an inflated cube placed inside a sphere and provides an almost uniform grid point distribution. The numerical schemes are based on a high-order variant of the Monotone Upstream-centered Schemes for Conservation Laws (MUSCL) pioneered by van Leer. In each cell the reconstructed left and right states are either obtained via a dimension-split piecewise-parabolic method or a piecewise-cubic reconstruction. The reconstructed states then serve as input to an approximate Riemann solver that determines the numerical fluxes at two Gaussian quadrature points along the cell boundary. The use of multiple quadrature points renders the resulting flux high-order. Three types of approximate Riemann solvers are compared, including the widely used solver of Rusanov, the solver of Roe and the new AUSM⁺-up solver of Liou that has been designed for low-Mach number flows. Spatial discretizations are paired with either a third-order or fourth-order total-variation-diminishing Runge–Kutta timestepping scheme to match the order of the spatial discretization. The numerical schemes are evaluated with several standard shallow-water test cases that emphasize accuracy and conservation properties. These tests show that the AUSM⁺-up flux provides the best overall accuracy, followed closely by the Roe solver. The Rusanov flux, with its simplicity, provides significantly larger errors by comparison. A brief discussion on extending the method to arbitrary order-of-accuracy is included.     

A general strategy for the optimization of Runge–Kutta schemes for wave propagation phenomena

We analyze optimized explicit Runge–Kutta schemes (RK) for computational aeroacoustics, and wave propagation phenomena in general. Exploiting the analysis developed in [S. Pirozzoli, Performance analysis and optimization of finite-difference schemes for wave propagation problems, J. Comput. Phys. 222 (2007) 809–831], we rigorously evaluate the performance of several time integration schemes in terms of appropriate error and cost metrics, and provide a general strategy to design Runge–Kutta methods tailored for specific applications. We present families of optimized second- and third-order Runge–Kutta schemes with up to seven stages, and describe their implementation in the framework of Williamson’s 2N$2N$-storage formulation [J.H. Williamson, Low-storage Runge–Kutta schemes, J. Comput. Phys. 35 (1980) 48–56]. Numerical simulations of the 1D linear advection equation and of the 2D linearized Euler equations are performed to demonstrate the validity of the theory and to quantify the improvement provided by optimized schemes.     

High order conservative Lagrangian schemes with Lax–Wendroff type time discretization for the compressible Euler equations

In this paper, we explore the Lax–Wendroff (LW) type time discretization as an alternative procedure to the high order Runge–Kutta time discretization adopted for the high order essentially non-oscillatory (ENO) Lagrangian schemes developed in 3 and 5. The LW time discretization is based on a Taylor expansion in time, coupled with a local Cauchy–Kowalewski procedure to utilize the partial differential equation (PDE) repeatedly to convert all time derivatives to spatial derivatives, and then to discretize these spatial derivatives based on high order ENO reconstruction. Extensive numerical examples are presented, for both the second-order spatial discretization using quadrilateral meshes [3] and third-order spatial discretization using curvilinear meshes [5]. Comparing with the Runge–Kutta time discretization procedure, an advantage of the LW time discretization is the apparent saving in computational cost and memory requirement, at least for the two-dimensional Euler equations that we have used in the numerical tests.     

Local DG method using WENO type limiters for convection–diffusion problems

The local discontinuous Galerkin (LDG) method is a spatial discretization procedure for convection–diffusion equations, which employs useful features from high resolution finite volume schemes, such as the exact or approximate Riemann solvers serving as numerical fluxes and limiters, which is termed as Runge–Kutta LDG (RKLDG) when TVD Runge–Kutta method is applied for time discretization. It has the advantage of flexibility in handling complicated geometry, h-p adaptivity, and efficiency of parallel implementation and has been used successfully in many applications. However, the limiters used to control spurious oscillations in the presence of strong shocks are less robust than the strategies of essentially non-oscillatory (ENO) and weighted ENO (WENO) finite volume and finite difference methods. In this paper, we investigated RKLDG methods with WENO and Hermite WENO (HWENO) limiters for solving convection–diffusion equations on unstructured meshes, with the goal of obtaining a robust and high order limiting procedure to simultaneously obtain uniform high order accuracy and sharp, non-oscillatory shock transition. Numerical results are provided to illustrate the behavior of these procedures.     

High-order parabolic beam approximation for aero-optics

The parabolic beam equations are solved using high-order compact differences for the Laplacians and Runge–Kutta integration along the beam path. The solution method is verified by comparison to analytical solutions for apertured beams and both constant and complex index of refraction. An adaptive 4th-order Runge–Kutta using an embedded 2nd-order method is presented that has demonstrated itself to be very robust. For apertured beams, the results show that the method fails to capture near aperture effects due to a violation of the paraxial approximation in that region. Initial results indicate that the problem appears to be correctable by successive approximations. A preliminary assessment of the effect of turbulent scales is undertaken using high-order Lagrangian interpolation. The results show that while high fidelity methods are necessary to accurately capture the large scale flow structure, the method may not require the same level of fidelity in sampling the density for the index of refraction. The solution is used to calculate a phase difference that is directly compared with that commonly calculated via the optical path difference. Propagation through a supersonic boundary layer shows that for longer wavelengths, the traditional method to calculate the optical path is less accurate than for shorter wavelengths. While unlikely to supplant more traditional methods for most aero-optics applications, the current method can be used to give a quantitative assessment of the other methods as well as being amenable to the addition of more physics.     

A high-order low-dispersion symmetry-preserving finite-volume method for compressible flow on curvilinear grids

A new high-order finite-volume method is presented that preserves the skew symmetry of convection for the compressible flow equations. The method is intended for Large-Eddy Simulations (LES) of compressible turbulent flows, in particular in the context of hybrid RANS–LES computations. The method is fourth-order accurate and has low numerical dissipation and dispersion. Due to the finite-volume approach, mass, momentum, and total energy are locally conserved. Furthermore, the skew-symmetry preservation implies that kinetic energy, sound-velocity, and internal energy are all locally conserved by convection as well. The method is unique in that all these properties hold on non-uniform, curvilinear, structured grids. Due to the conservation of kinetic energy, there is no spurious production or dissipation of kinetic energy stemming from the discretization of convection. This enhances the numerical stability and reduces the possible interference of numerical errors with the subgrid-scale model. By minimizing the numerical dispersion, the numerical errors are reduced by an order of magnitude compared to a standard fourth-order finite-volume method.     

Adaptive hierarchic transformations for dynamically p-enriched slope-limiting over discontinuous Galerkin systems of generalized equations

We study a family of generalized slope limiters in two dimensions for Runge–Kutta discontinuous Galerkin (RKDG) solutions of advection-diffusion systems. We analyze the numerical behavior of these limiters applied to a pair of model problems, comparing the error of the approximate solutions, and discuss each limiter’s advantages and disadvantages. We then introduce a series of coupled p-enrichment schemes that may be used as standalone dynamic p-enrichment strategies, or may be augmented via any in the family of variable-in-p slope limiters presented.     

An adaptive implicit–explicit scheme for the DNS and LES of compressible flows on unstructured grids

An adaptive implicit–explicit scheme for Direct Numerical Simulation (DNS) and Large-Eddy Simulation (LES) of compressible turbulent flows on unstructured grids is developed. The method uses a node-based finite-volume discretization with Summation-by-Parts (SBP) property, which, in conjunction with Simultaneous Approximation Terms (SAT) for imposing boundary conditions, leads to a linearly stable semi-discrete scheme. The solution is marched in time using an Implicit–Explicit Runge–Kutta (IMEX-RK) time-advancement scheme. A novel adaptive algorithm for splitting the system into implicit and explicit sets is developed. The method is validated using several canonical laminar and turbulent flows. Load balance for the new scheme is achieved by a dual-constraint, domain decomposition algorithm. The scalability and computational efficiency of the method is investigated, and memory savings compared with a fully implicit method is demonstrated. A notable reduction of computational costs compared to both fully implicit and fully explicit schemes is observed.     

An interface treating technique for compressible multi-medium flow with Runge–Kutta discontinuous Galerkin method

The high-order accurate Runge–Kutta discontinuous Galerkin (RKDG) method is applied to the simulation of compressible multi-medium flow, generalizing the interface treating method given in Chertock et al. (2008) [9]. In mixed cells, where the interface is located, Riemann problems are solved to define the states on both sides of the interface. The input states to the Riemann problem are obtained by extrapolation to the cell boundary from solution polynomials in the neighbors of the mixed cell. The level set equation is solved by using a high-order accurate RKDG method for Hamilton–Jacobi equations, resulting in a unified DG solver for the coupled problem. The method is conservative if we include the states in the mixed cells, which are however not used in the updating of the numerical solution in other cells. The states in the mixed cells are plotted to better evaluate the conservation errors, manifested by overshoots/undershoots when compared with states in neighboring cells. These overshoots/undershoots in mixed cells are problem dependent and change with time. Numerical examples show that the results of our scheme compare well with other methods for one and two-dimensional problems. In particular, the algorithm can capture well complex flow features of the one-dimensional shock entropy wave interaction problem and two-dimensional shock–bubble interaction problem.     

Higher Order KFVS Algorithms Using Compact Upwind Difference Operators

A family of high order accurate compact upwind difference operators have been used, together with the split fluxes of the KFVS (kinetic flux vector splitting) scheme to obtain high order semidiscretizations of the 2D Euler equations of inviscid gas dynamics in general coordinates. A TVD multistage Runge–Kutta time stepping scheme is used to compute steady states for selected transonic/supersonic flow problems which indicate the higher accuracy and low diffusion realizable in such schemes.     

A high order kinetic flux-vector splitting method for the reduced five-equation model of compressible two-fluid flows

We present a high order kinetic flux-vector splitting (KFVS) scheme for the numerical solution of a conservative interface-capturing five-equation model of compressible two-fluid flows. This model was initially introduced by Wackers and Koren (2004) [21]. The flow equations are the bulk equations, combined with mass and energy equations for one of the two fluids. The latter equation contains a source term in order to account for the energy exchange. We numerically investigate both one- and two-dimensional flow models. The proposed numerical scheme is based on the direct splitting of macroscopic flux functions of the system of equations. In two space dimensions the scheme is derived in a usual dimensionally split manner. The second order accuracy of the scheme is achieved by using MUSCL-type initial reconstruction and Runge–Kutta time stepping method. For validation, the results of our scheme are compared with those from the high resolution central scheme of Nessyahu and Tadmor [14]. The accuracy, efficiency and simplicity of the KFVS scheme demonstrate its potential for modeling two-phase flows.     

Hybrid weighted essentially non-oscillatory schemes with different indicators

A key idea in finite difference weighted essentially non-oscillatory (WENO) schemes is a combination of lower order fluxes to obtain a higher order approximation. The choice of the weight to each candidate stencil, which is a nonlinear function of the grid values, is crucial to the success of WENO schemes. For the system case, WENO schemes are based on local characteristic decompositions and flux splitting to avoid spurious oscillation. But the cost of computation of nonlinear weights and local characteristic decompositions is very high. In this paper, we investigate hybrid schemes of WENO schemes with high order up-wind linear schemes using different discontinuity indicators and explore the possibility in avoiding the local characteristic decompositions and the nonlinear weights for part of the procedure, hence reducing the cost but still maintaining non-oscillatory properties for problems with strong shocks. The idea is to identify discontinuity by an discontinuity indicator, then reconstruct numerical flux by WENO approximation in discontinuous regions and up-wind linear approximation in smooth regions. These indicators are mainly based on the troubled-cell indicators for discontinuous Galerkin (DG) method which are listed in the paper by Qiu and Shu (J. Qiu, C.-W. Shu, A comparison of troubled-cell indicators for Runge–Kutta discontinuous Galerkin methods using weighted essentially non-oscillatory limiters, SIAM Journal of Scientific Computing 27 (2005) 995–1013). The emphasis of the paper is on comparison of the performance of hybrid scheme using different indicators, with an objective of obtaining efficient and reliable indicators to obtain better performance of hybrid scheme to save computational cost. Detail numerical studies in one- and two-dimensional cases are performed, addressing the issues of efficiency (less CPU time and more accurate numerical solution), non-oscillatory property.