High-order Central Difference Research Articles

On the use of high order central difference schemes for differential equation based wall distance computations

A computationally efficient high-order solver is developed to compute the wall distances by solving the relevant partial differential equations, namely: Eikonal, Hamilton–Jacobi (HJ) and Poisson equations. In contrast to the upwind schemes widely used in the literature, we explore the suitability of high-order central difference schemes (explicit/compact) for the wall-distance computation. While solving the Hamilton–Jacobi equation, the high-order central difference schemes performed approximately 1.4–2.8 times faster than the upwind schemes with a marginal improvement in the solution accuracy. A new pseudo HJ formulation based on the localized artificial diffusivity (LAD) approach has been proposed. It is demonstrated to predict results with an accuracy comparable to that of the Eikonal equation and the simulations are ≈ 1.5 times faster than the baseline HJ solver using upwind schemes. A curvature correction is also incorporated in the HJ equation to correct for the near-wall errors due to concave/convex wall curvatures. We demonstrate the efficacy of the proposed methods on both the steady and unsteady test cases and exploit the unsteady wall-distance solver to estimate the instantaneous shape and burning surface area of a dendrite propellant grain in a solid propellant rocket motor.

Numerical visualisation on applying a five-stage multi-scaled square turbulence-generating grid to generating turbulence

This study presents the numerical visualisation of turbulence generated by a five-stage square-type multi-scale turbulence grid. The present study applies a turbulence grid with an increased number of stages to five, in contrast to previous studies using multi-scale turbulence grids with four stages. For the present numerical visualisation, direct numerical analysis is used in this analysis. This analysis uses the high-order central difference schemes and the third-order Runge-Kutta scheme. In addition to the multi-scale turbulence grid with five stages, the present study also covers two turbulence grids with four stages for reference. By using the external force term in the governing equations, the present multi-scale turbulence grids were formed numerically in the computational domain. The spatial homogeneity of the downstream turbulence field is approached in this work. The turbulence generated by the turbulence grids with four stages is characterised by a flow structure with a doughnut-like instantaneous velocity figure in the downstream region. By examining contour diagrams of the streamwise instantaneous velocity, the present study can see that the spatial homogeneity of the turbulence field downstream is improved by increasing the number of stages for the turbulence grid.

A FFT accelerated high order finite difference method for elliptic boundary value problems over irregular domains

For elliptic boundary value problems (BVPs) involving irregular domains and Robin boundary condition, no numerical method is known to deliver a fourth order convergence and O(Nlog⁡N) efficiency, where N stands for the total degree-of-freedom of the system. Based on the matched interface and boundary (MIB) and fast Fourier transform (FFT) schemes, a new finite difference method is introduced for such problems, which involves two main components. First, a ray-casting MIB scheme is proposed to handle different types of boundary conditions, including Dirichlet, Neumann, Robin, and their mix combinations. By enclosing the concerned irregular domain by a large enough cubic domain, the ray-casting MIB scheme generates necessary fictitious values outside the irregular domain by imposing boundary conditions along the normal direction of the boundary, so that a high order central difference discretization of the Laplacian can be formed. Second, an augmented MIB formulation is built, in which Cartesian derivative jumps are reconstructed on the boundary as auxiliary variables. By treating such variables as unknowns, the discrete Laplacian can be efficiently inverted by the FFT algorithm, in the Schur complement solution of the augmented system. The accuracy and efficiency of the proposed augmented MIB method are numerically examined by considering various elliptic BVPs in two and three dimensions. Numerical results indicate that the new algorithm not only achieves a fourth order of accuracy in treating irregular domains and complex boundary conditions, but also maintains the FFT efficiency.

FFT-Based High Order Central Difference Schemes for Poisson’s Equation with Staggered Boundaries

This work concerns with the development of fast and high order algorithms for solving a single variable Poisson’s equation with rectangular domains and uniform meshes, but involving staggered boundaries. Here the staggered boundary means that the boundary is located midway between two adjacent grid nodes. Due to the popularity of staggered grids in scientific computing for solving multiple variables partial differential equations (PDEs), the planned development deserves further studies, but is rarely reported in the literature, because grand challenges exist for spectral methods, compact finite differences, and fast Fourier transform (FFT) algorithms in handling staggered boundaries. A systematic approach is introduced in this paper to attack various open problems in this regard, which is a natural generalization of a recently developed Augmented Matched Interface and Boundary (AMIB) method for non-staggered boundaries. Formulated through immersed boundary problems with zero-padding solutions, the AMIB method combines arbitrarily high order central differences with the FFT inversion. Over staggered boundaries, the proposed AMIB method can handle Dirichlet, Neumann, Robin or any combination of boundary conditions. Convergence orders in four, six and eight are numerically validated for the AMIB method in both two and three dimensions. Moreover, the proposed AMIB method performs well for some challenging problems, such as low regularity solution near boundary, PDE solution not satisfying the boundary condition, and involving both staggered and non-staggered boundaries on two ends.

FFT-based high order central difference schemes for three-dimensional Poisson's equation with various types of boundary conditions

In this paper, a unified approach is introduced to implement high order central difference schemes for solving Poisson's equation via the fast Fourier transform (FFT). Popular high order fast Poisson solvers in the literature include compact finite differences and spectral methods. However, FFT-based high order central difference schemes have never been developed for Poisson problems, because with long stencils, central differences require fictitious nodes outside the boundary, which poses a challenge to integrate boundary conditions in FFT computations. To overcome this difficulty, several layers of exterior grid lines are introduced to convert the problem to an immersed boundary problem with zero-padding solutions beyond the original cubic domain. Over the boundary of the enlarged cubic domain, the anti-symmetric property is naturally satisfied so that the FFT fast inversion is feasible, while the immersed boundary problem can be efficiently solved by the proposed augmented matched interface and boundary (AMIB) method. As the first fast Poisson solver based on high order central differences, the AMIB method can be easily implemented in any dimension, due to its tensor product nature of the discretization. As a systematical approach, the AMIB method can be made to arbitrarily high order in principle, and can handle the Dirichlet, Neumann, Robin or any combination of boundary conditions. The accuracy, efficiency, and robustness of the proposed AMIB method are numerically validated by considering various Poisson problems in two and three dimensions.

High-order accurate large-eddy simulations of compressible viscous flow in cylindrical coordinates

Dynamic interactions of shock waves and turbulent structures occur in a wide range of applications. The accurate simulation of turbulence in these flows requires a numerical scheme with minimal dissipation whereas the shock-capturing requires a dissipative scheme to stabilise the solution in the shock-wave discontinuous regions. These contradictory requirements make simulations of these flows extremely challenging. A compatible implementation of the solid boundary conditions at the point of incepting the acoustic waves and their internalisation into hydrodynamic instabilities is also challenging.Here we present, a high-fidelity, massively parallel and accurate solver for both direct numerical simulations and large-eddy simulations of supersonic turbulent flows in cylindrical coordinates. A hybrid WENO/ high order central difference scheme is used for the spatial discretisation with a fourth-order five-stage 2N-storage Runge–Kutta for the time integration. The least square contraction of Lilly [1] is utilised for large-eddy subgrid-scale modelling. The solid surface boundary condition is implemented using the offset wall and ghost cell technique.The numerical implementation is validated through a wide range of test cases from simple to more complex cases. The numerical results show good agreement with analytical solutions and available experimental results. In the large-eddy simulation of a supersonic under-expanded impinging jet, it is observed that the conventional Ducros sensor leads to a more dissipative hybrid scheme. In this paper, a new shock sensor based on a WENO smoothness indicator is proposed and it is demonstrated to perform better especially in the simulation of the supersonic under-expanded impinging jet by limiting the expensive WENO scheme to the discontinuous regions and reducing the computational cost by approximately 10%.

Entropy Stable Method for the Euler Equations Revisited: Central Differencing via Entropy Splitting and SBP

The two decades old high order central differencing via entropy splitting and summation-by-parts (SBP) difference boundary closure of Olsson and Oliger (Energy and maximum norm estimates for nonlinear conservation laws, 1994), Gerritsen and Olsson (J Comput Phys 129:245–262, 1996) and Yee et al. (J Comput Phys 162:33–81, 2000) is revisited. The objective of this paper is to prove for the first time that the entropy split methods based on physical entropies are entropy stable methods for central differencing with SBP operators for both periodic and non-periodic boundary conditions for nonlinear Euler equations. Standard high order spatial central differencing as well as high order central spatial dispersion relation preserving spatial differencing is part of the entropy stable methodology framework. The proof is to replace the spatial derivatives by SBP difference operators in the entropy split form of the equations using the physical entropy of the Euler equations. The numerical boundary closure follows directly from the SBP operator. No additional numerical boundary procedure is required. In contrast, Tadmor-type entropy conserving methods (Acta Numer 12:451–512, 2003) using mathematical entropies and more recently the methods in Winters and Gassner (J Comput Phys 304:72–108, 2016), do not naturally come with a numerical boundary closure and a generalized SBP operator has to be developed (Roanocha in Generalized summation-by-parts operators and variable coefficients, 2018. arXiv:1705.10541v2 [math.NA]). Long time integration of 2D and 3D test cases is included to show the comparison of this efficient entropy stable method with the Tadmor-type of entropy conservative methods. Studies also include the comparison among the three skew-symmetric splittings on their nonlinear stability and accuracy performance without added numerical dissipations for smooth flows. These are, namely, entropy splitting, Ducros et al. splitting and the Kennedy and Grubber splitting.

Assessment of artificial fluid properties for high-order accurate large-eddy simulations of shock-free compressible turbulent flows with strong temperature gradients

The capability of artificial fluid properties (AFP) to serve as a sub-grid scale (SGS) model and to provide stability in high-order accurate large-eddy simulations (LES) of shock-free compressible turbulent flows with strong temperature gradients is assessed. Predictive LES of such flows, which are directly relevant to applications such as engines, require an accurate description of turbulence, which typically makes them susceptible to issues of accuracy and numerical instabilities. In the present implementation, AFP schemes along with a high-order compact filter ensure robustness of otherwise non-dissipative high-order central difference spatial discretisation schemes. It is shown that a judicious choice of both filter and AFP parameters is required in order for AFP to be used as an accurate SGS model. Optimal parameter values for the AFP and the corresponding filter parameter are identified and validated in the framework of the current LES solver. First, an inviscid Taylor-Green vortex configuration is examined to evaluate the ideal combination of the filter parameter and the coefficients in the AFP framework. It is verified that AFP works as an alternative model for the unclosed subgrid-scale (SGS) dissipation terms. The AFP results are also compared to those from the dynamic Smagorinsky model (DSM), with good agreement. Second, the optimal parameter values of AFP obtained from the inviscid Taylor-Green vortex tests are validated on a series of test configurations with increasing complexity, that include an inviscid 1D advection problem involving a thermal contact discontinuity, an inviscid 2D Kelvin-Helmholtz instability test and 3D experimental unheated and heated spatially developing turbulent round jets. Robust and accurate predictions are shown with the LES solver employing the optimal model parameters in each of these tests. In particular, good agreement between the LES predictions and experimental measurements for both the flow and temperature statistics is achieved for the laboratory jets, demonstrating the accuracy of the LES/AFP simulations applied to practical turbulent configurations. The robustness of the solver in the presence of strong temperature gradients is also demonstrated by simulating a round jet case having very high temperature gradients similar to those occurring in reacting flows. In contrast, the LES/DSM simulation of this case is noted to be unstable, and failed to provide converged results.

A high-order accurate AUSM$$^+$$-up approach for simulations of compressible multiphase flows with linear viscoelasticity

An Eulerian approach for simulations of wave propagation in multiphase, viscoelastic media is developed in the context of the Advection Upstream Splitting Method (AUSM). We extend the AUSM scheme to the five-equation model for simulations of interfaces between gases, liquids, and solids with constitutive relations appropriately transported. In this framework, the solid’s deformations are assumed to be infinitesimally small such that they can be modeled using linear viscoelastic models, e.g., generalized Zener. The Eulerian framework addresses the challenge of calculating strains, more naturally expressed in a Lagrangian framework, by using a hypoelastic model that takes an objective Lie derivative of the constitutive relation to transform strains into velocity gradients. Our approach introduces elastic stresses in the convective fluxes that are treated by generalizing AUSM flux-vector splitting (FVS) to account for the Cauchy stress tensor. We determine an appropriate discretization of non-conservative equations that appear in the five-equation multiphase model with AUSM schemes to prevent spurious oscillations at material interfaces. The framework’s spatial scheme is solution adaptive with a discontinuity sensor discriminating between smooth and discontinuous regions. The smooth regions are computed using explicit high-order central differences. At discontinuous regions (i.e., shocks, material interfaces, and contact surfaces), the convective fluxes are treated using a high-order Weighted Essentially Non-Oscillatory (WENO) scheme with $$\hbox {AUSM}^+$$ -up for upwinding. The framework is used to simulate one-dimensional (1D) and two-dimensional (2D) problems that demonstrate the ability to maintain equilibrium interfacial conditions and solve challenging multi-dimensional and multi-material problems.

A high-order accurate five-equations compressible multiphase approach for viscoelastic fluids and solids with relaxation and elasticity

A novel Eulerian approach is proposed for numerical simulations of wave propagation in viscoelastic media, for application to shocks interacting with interfaces between fluids and solids. We extend the five-equations multiphase interface-capturing model, based on the idea that all the materials (gases, liquids, solids) obey the same equation of state with spatially varying properties, to incorporate the desired constitutive relation; in this context, interfaces are represented by discontinuities in material properties. We consider problems in which the deformations are small, such that the substances can be described by linear constitutive relations, specifically, Maxwell, Kelvin–Voigt or generalized Zener models. The main challenge lies in representing the combination of viscoelastic, multiphase and compressible flow. One particular difficulty is the calculation of strains in an Eulerian framework, which we address by using a conventional hypoelastic model in which an objective time derivative (Lie derivative) of the constitutive relation is taken to evolve strain rates instead. The resulting eigensystem is analyzed to identify wave speeds and characteristic variables. The spatial scheme is based on a solution-adaptive formulation, in which a discontinuity sensor discriminates between smooth and discontinuous regions. To compute the convective fluxes, explicit high-order central differences are applied in smooth regions, while a high-order finite-difference Weighted Essentially Non-Oscillatory (WENO) scheme is used at discontinuities (shocks, material interfaces and contacts). The numerical method is verified in a comprehensive fashion using a series of smooth and discontinuous (shocks and interfaces), one- and two-dimensional test problems.

High-order central difference scheme for Caputo fractional derivative

In this paper we propose a class of central difference schemes for resolving the Caputo fractional derivative. The accuracy may reach any selected integer order. More precisely, the Caputo fractional derivative operator is decomposed into symmetric and antisymmetric components. Starting from difference schemes of lower order accuracy for each component, we enhance the accuracy by a weighted average of shifted differences. The weights are calculated by matching the symbols of the scheme and the operators. We further illustrate the application of the proposed schemes to a fractional advection–diffusion equation. Together with the Crank–Nicolson algorithm, it reaches designed accuracy order, and is unconditionally stable. Numerical tests are presented to demonstrate the nice features.

Stability Analysis of the Inverse Lax–Wendroff Boundary Treatment for High Order Central Difference Schemes for Diffusion Equations

In this paper, high order central finite difference schemes in a finite interval are analyzed for the diffusion equation. Boundary conditions of the initial-boundary value problem are treated by the simplified inverse Lax–Wendroff procedure. For the fully discrete case, a third order explicit Runge–Kutta method is used as an example for the analysis. Stability is analyzed by both the Gustafsson, Kreiss and Sundstrom theory and the eigenvalue visualization method on both semi-discrete and fully discrete schemes. The two different analysis techniques yield consistent results. Numerical tests are performed to demonstrate and validate the analysis results.

Application of a Multi-dimensional Limiting Process to Central-Upwind Schemes for Solving Hyperbolic Systems of Conservation Laws

In this paper, we study semi-discrete central-upwind difference schemes with a modified multi-dimensional limiting process (MLP) to solve two-dimensional hyperbolic systems of conservation laws. In general, high-order central difference schemes for conservation laws involve no Riemann solvers or characteristic decompositions but have a tendency to smear linear discontinuities. To overcome this drawback of central-upwind schemes, we use a MLP that uses multi-dimensional information for slope limitation to control the oscillations across discontinuities for multi-dimensional applications. Some numerical results are provided to demonstrate the performance of the proposed scheme.

Wave-equation-based travel-time seismic tomography – Part 1: Method

Abstract. In this paper, we propose a wave-equation-based travel-time seismic tomography method with a detailed description of its step-by-step process. First, a linear relationship between the travel-time residual Δt = Tobs–Tsyn and the relative velocity perturbation δ c(x)/c(x) connected by a finite-frequency travel-time sensitivity kernel K(x) is theoretically derived using the adjoint method. To accurately calculate the travel-time residual Δt, two automatic arrival-time picking techniques including the envelop energy ratio method and the combined ray and cross-correlation method are then developed to compute the arrival times Tsyn for synthetic seismograms. The arrival times Tobs of observed seismograms are usually determined by manual hand picking in real applications. Travel-time sensitivity kernel K(x) is constructed by convolving a~forward wavefield u(t,x) with an adjoint wavefield q(t,x). The calculations of synthetic seismograms and sensitivity kernels rely on forward modeling. To make it computationally feasible for tomographic problems involving a large number of seismic records, the forward problem is solved in the two-dimensional (2-D) vertical plane passing through the source and the receiver by a high-order central difference method. The final model is parameterized on 3-D regular grid (inversion) nodes with variable spacings, while model values on each 2-D forward modeling node are linearly interpolated by the values at its eight surrounding 3-D inversion grid nodes. Finally, the tomographic inverse problem is formulated as a regularized optimization problem, which can be iteratively solved by either the LSQR solver or a~nonlinear conjugate-gradient method. To provide some insights into future 3-D tomographic inversions, Fréchet kernels for different seismic phases are also demonstrated in this study.

Consistent numerical diffusion terms for simulating compressible multicomponent flows

An interface-capturing method using a high-order central difference scheme is presented for simulations of compressible multicomponent flows. The present method adds consistent numerical diffusion terms to robustly capture interface discontinuities, while maintaining the velocity, pressure, and temperature equilibriums at interfaces. Analysis of the numerical errors generated at the interfaces leads to a proposal of new consistent numerical diffusion terms. The method solves a fully conservative form of the total mass, momentum, total energy, and species-mass, and an additional advection form of the specific heats ratio to preserve the pressure oscillation-free property at the interfaces. Several one- and two-dimensional problems are used to verify the proposed method.

Stability of Finite Difference Discretizations of Multi-Physics Interface Conditions

AbstractWe consider multi-physics computations where the Navier-Stokes equations of compressible fluid flow on some parts of the computational domain are coupled to the equations of elasticity on other parts of the computational domain. The different subdomains are separated by well-defined interfaces. We consider time accurate computations resolving all time scales. For such computations, explicit time stepping is very efficient. We address the issue of discrete interface conditions between the two domains of different physics that do not lead to instability, or to a significant reduction of the stable time step size. Finding such interface conditions is non-trivial.We discretize the problem with high order centered difference approximations with summation by parts boundary closure. We derive L2 stable interface conditions for the linearized one dimensional discretized problem. Furthermore, we generalize the interface conditions to the full non-linear equations and numerically demonstrate their stable and accurate performance on a simple model problem. The energy stable interface conditions derived here through symmetrization of the equations contain the interface conditions derived through normal mode analysis by Banks and Sjögreen in [8] as a special case.

Stability of molecular dynamics simulations of classical systems

The existence of a shadow Hamiltonian H̃ for discrete classical dynamics, obtained by an asymptotic expansion for a discrete symplectic algorithm, is employed to determine the limit of stability for molecular dynamics (MD) simulations with respect to the time-increment h of the discrete dynamics. The investigation is based on the stability of the shadow energy, obtained by including the first term in the asymptotic expansion, and on the exact solution of discrete dynamics for a single harmonic mode. The exact solution of discrete dynamics for a harmonic potential with frequency ω gives a criterion for the limit of stability h ≤ 2/ω. Simulations of the Lennard-Jones system and the viscous Kob-Andersen system show that one can use the limit of stability of the shadow energy or the stability criterion for a harmonic mode on the spectrum of instantaneous frequencies to determine the limit of stability of MD. The method is also used to investigate higher-order central difference algorithms, which are symplectic and also have shadow Hamiltonians, and for which one can also determine the exact criteria for the limit of stability of a single harmonic mode. A fourth-order central difference algorithm gives an improved stability with a factor of 3, but the overhead of computer time is a factor of at least two. The conclusion is that the second-order "Verlet"-algorithm, most commonly used in MD, is superior. It gives the exact dynamics within the limit of the asymptotic expansion and this limit can be estimated either from the conserved shadow energy or from the instantaneous spectrum of harmonic modes.

On the Use of Temporal Reachout Technique for Characteristics Method with Time-Line Cubic Spline Interpolation

ABSTRACTThe study had indicated that the computational performances of the characteristics method with the time-line cubic spline interpolation are related to the endpoint constraint, especially for large Courant number in which the foot of the characteristic trajectory is located near the endpoint. The first derivative endpoint constraint with higher-order central difference approximation provides better simulation results among various endpoint constraints, but it still induces some degree of numerical error. In this study, by locating the foot of the characteristic trajectory away from the endpoint, the temporal reachout technique is proposed to avoid the effect of endpoint constraint on the time-line cubic spline interpolation. Modeling the transport of a Gaussian concentration distribution in a uniform flow with constant diffusion coefficient and the viscous Burgers equation is used to examine the temporal reachout technique. The outcomes show that the temporal reachout technique yields much better simulation results than the first derivative endpoint constraint with higher-order central difference approximation. The effect of endpoint constraint on the time-line cubic spline interpolation can be greatly diminished by the use of the temporal reachout technique.

THE PARAMETER-ROBUST NUMERICAL METHOD BASED ON DEFECT-CORRECTION TECHNIQUE FOR SINGULARLY PERTURBED DELAY DIFFERENTIAL EQUATIONS WITH LAYER BEHAVIOR

In this article, we consider a defect-correction method based on finite difference scheme for solving a singularly perturbed delay differential equation. We solve the equation using upwind finite difference scheme on piecewise-uniform Shishkin mesh, then apply the defect-correction technique that combines the stability of the upwind scheme and the higher-order central difference scheme. The method is shown to be convergent uniformly in the perturbation parameter and almost second-order convergence measured in the discrete maximum norm is obtained. Numerical results are presented, which are in agreement with the theoretical findings.

Construction of low dissipative high-order well-balanced filter schemes for non-equilibrium flows

The goal of this paper is to generalize the well-balanced approach for non-equilibrium flow studied by Wang et al. (2009) [29] to a class of low dissipative high-order shock-capturing filter schemes and to explore more advantages of well-balanced schemes in reacting flows. More general 1D and 2D reacting flow models and new examples of shock turbulence interactions are provided to demonstrate the advantage of well-balanced schemes. The class of filter schemes developed by Yee et al. (1999) [33], Sjögreen and Yee (2004) [27] and Yee and Sjögreen (2007) [38] consist of two steps, a full time step of spatially high-order non-dissipative base scheme and an adaptive non-linear filter containing shock-capturing dissipation. A good property of the filter scheme is that the base scheme and the filter are stand-alone modules in designing. Therefore, the idea of designing a well-balanced filter scheme is straightforward, i.e. choosing a well-balanced base scheme with a well-balanced filter (both with high-order accuracy). A typical class of these schemes shown in this paper is the high-order central difference schemes/predictor–corrector (PC) schemes with a high-order well-balanced WENO filter. The new filter scheme with the well-balanced property will gather the features of both filter methods and well-balanced properties: it can preserve certain steady-state solutions exactly; it is able to capture small perturbations, e.g. turbulence fluctuations; and it adaptively controls numerical dissipation. Thus it shows high accuracy, efficiency and stability in shock/turbulence interactions. Numerical examples containing 1D and 2D smooth problems, 1D stationary contact discontinuity problem and 1D turbulence/shock interactions are included to verify the improved accuracy, in addition to the well-balanced behavior.