High-order Accurate Finite Difference Research Articles

Compatible diagonal-norm staggered and upwind SBP operators

The main motivation with the present study is to achieve a provably stable high-order accurate finite difference discretisation of linear first-order hyperbolic problems on a staggered grid. The use of a staggered grid makes it non-trivial to discretise advective terms. To overcome this difficulty we discretise the advective terms using upwind Summation-By-Parts (SBP) operators, while the remaining terms are discretised using staggered SBP operators. The upwind and staggered SBP operators (for each order of accuracy) are compatible, here meaning that they are based on the same diagonal norms, allowing for energy estimates to be formulated. The boundary conditions are imposed using a penalty (SAT) technique, to guarantee linear stability. The resulting SBP-SAT approximations lead to fully explicit ODE systems. The accuracy and stability properties are demonstrated for linear hyperbolic problems in 1D, and for the 2D linearised Euler equations with constant background flow. The newly derived upwind and staggered SBP operators lead to significantly more accurate numerical approximations, compared with the exclusive usage of (previously derived) central-difference first derivative SBP operators.

High order accurate finite difference schemes based on symmetry preservation

In this paper, we present a mathematical approach that is based on modified equations and the method of equivariant moving frames for construction of high order accurate invariant finite difference schemes that preserve Lie symmetry groups of underlying partial differential equations (PDEs). In the proposed approach, invariant (or symmetry preserving) numerical schemes with a desired (or fixed) order of accuracy are constructed from some non-invariant (base) numerical schemes. Modified forms of PDEs are used to improve the order of accuracy of existing schemes and these modified forms are obtained through addition of defect correction terms to the original forms of PDEs. These defect correction terms of modified PDEs that are noted from truncation error analysis are either completely removed from schemes or their representation is significantly simplified by considering convenient moving frames. This feature of the proposed method can especially be useful to avoid cumbersome numerical representations when high order schemes are developed from low order ones via the method of modified equations. The proposed method is demonstrated via construction of invariant numerical schemes with fixed (and higher) order of accuracy for some common linear and nonlinear problems (including the linear advection–diffusion equation in 1D and 2D, inviscid Burgers' equation, and viscous Burgers' equation) and the performance of these invariant numerical schemes is further evaluated. Our results indicate that such invariant numerical schemes obtained from existing base numerical schemes have the potential to significantly improve the quality of results not only in terms of desired higher order accuracy but also in the context of preservation of appropriate symmetry properties of underlying PDEs.

High Order Accurate Finite Difference Modeling of Seismo-Acoustic Wave Propagation in a Moving Atmosphere and a Heterogeneous Earth Model Coupled Across a Realistic Topography

We develop a numerical method for simultaneously simulating acoustic waves in a realistic moving atmosphere and seismic waves in a heterogeneous earth model, where the motions are coupled across a realistic topography. We model acoustic wave propagation by solving the linearized Euler equations of compressible fluid mechanics. The seismic waves are modeled by the elastic wave equation in a heterogeneous anisotropic material. The motion is coupled by imposing continuity of normal velocity and normal stresses across the topographic interface. Realistic topography is resolved on a curvilinear grid that follows the interface. The governing equations are discretized using high order accurate finite difference methods that satisfy the principle of summation by parts. We apply the energy method to derive the discrete interface conditions and to show that the coupled discretization is stable. The implementation is verified by numerical experiments, and we demonstrate a simulation of coupled wave propagation in a windy atmosphere and a realistic earth model with non-planar topography.

A Family of Compact Finite Difference Formulations for Three-Space Dimensional Nonlinear Poisson’s Equations in Cartesian Coordinates

In this article, a family of high order accurate compact finite difference scheme for obtaining the approximate solution values of mildly nonlinear elliptic boundary value problems in three-space dimensions has been developed. The discretization formula is developed on a non-uniform meshes, which helps in resolving boundary and/or interior layers. The scheme involves 27 points single computational cell to achieve high order truncation errors. The proposed scheme has been applied to solve convection–diffusion equation, Helmholtz equation and nonlinear Poisson’s equation. A detailed convergence theory for the new compact scheme has been proposed using irreducible and monotone property of the iteration matrix. Numerical results show that the new compact scheme exhibit better performance in terms of \(l^\infty \)- and \(l^2\)-error of the exact and approximate solution values.

On the reduction of aerofoil–turbulence interaction noise associated with wavy leading edges

An aerofoil leading-edge profile based on wavy (sinusoidal) protuberances/tubercles is investigated to understand the mechanisms by which they are able to reduce the noise produced through the interaction with turbulent mean flow. Numerical simulations are performed for non-lifting flat-plate aerofoils with straight and wavy leading edges (denoted by SLE and WLE, respectively) subjected to impinging turbulence that is synthetically generated in the upstream zone (free-stream Mach number of 0.24). Full three-dimensional Euler (inviscid) solutions are computed for this study thereby eliminating self-noise components. A high-order accurate finite-difference method and artefact-free boundary conditions are used in the current simulations. Various statistical analysis methods, including frequency spectra, are implemented to aid the understanding of the noise-reduction mechanisms. It is found with WLEs, unlike the SLE, that the surface pressure fluctuations along the leading edge exhibit a significant source-cutoff effect due to geometric obliqueness which leads to reduced levels of radiated sound pressure. It is also found that there exists a phase interference effect particularly prevalent between the peak and the hill centre of the WLE geometry, which contributes to the noise reduction in the mid- to high-frequency range.

Parametric numerical study of electrokinetic instability in cross-shaped microchannels

A multiblock in-house solver based on high-order accurate finite-difference schemes is developed to simulate electrokinetic instability (EKI) phenomena in realistic cross-shaped microchannels. Parametric numerical studies are systematically carried out to investigate the effects of different parameters on EKI including the conductivity ratio, the electric field strength, the channel depth, and the electric field ratios. Simulation results show that the solver has the capability to capture the threshold electric field for the onset of instabilities and it offers a good prediction for the critical features of EKI compared to the experiments presented in the literature. In addition, from the general parametric studies several useful guidelines are explored. Such guidelines can be applied to enhance the mixing efficiency in electrokinetic microdevices.

High-fidelity numerical simulation of solitons in the nerve axon

High-order accurate finite difference schemes are derived for a non-linear soliton model of nerve signal propagation in axons. Two types of well-posed boundary conditions are analysed. The boundary closures are based on the summation-by-parts (SBP) framework and the boundary conditions are imposed using a penalty (SAT) technique, to guarantee linear stability. The resulting SBP-SAT approximation is time-integrated with an explicit finite difference method. The accuracy and stability properties of the newly derived finite difference approximations are demonstrated for an analytic soliton solution.

Dynamic earthquake rupture simulations on nonplanar faults embedded in 3D geometrically complex, heterogeneous elastic solids

Dynamic propagation of shear ruptures on a frictional interface in an elastic solid is a useful idealization of natural earthquakes. The conditions relating discontinuities in particle velocities across fault zones and tractions acting on the fault are often expressed as nonlinear friction laws. The corresponding initial boundary value problems are both numerically and computationally challenging. In addition, seismic waves generated by earthquake ruptures must be propagated for many wavelengths away from the fault. Therefore, reliable and efficient numerical simulations require both provably stable and high order accurate numerical methods.We present a high order accurate finite difference method for: a) enforcing nonlinear friction laws, in a consistent and provably stable manner, suitable for efficient explicit time integration; b) dynamic propagation of earthquake ruptures along nonplanar faults; and c) accurate propagation of seismic waves in heterogeneous media with free surface topography.We solve the first order form of the 3D elastic wave equation on a boundary-conforming curvilinear mesh, in terms of particle velocities and stresses that are collocated in space and time, using summation-by-parts (SBP) finite difference operators in space. Boundary and interface conditions are imposed weakly using penalties. By deriving semi-discrete energy estimates analogous to the continuous energy estimates we prove numerical stability. The finite difference stencils used in this paper are sixth order accurate in the interior and third order accurate close to the boundaries. However, the method is applicable to any spatial operator with a diagonal norm satisfying the SBP property. Time stepping is performed with a 4th order accurate explicit low storage Runge-Kutta scheme, thus yielding a globally fourth order accurate method in both space and time. We show numerical simulations on band limited self-similar fractal faults revealing the complexity of rupture dynamics on rough faults.

Boundary conditions and stability of a perfectly matched layer for the elastic wave equation in first order form

In computations, it is now common to surround artificial boundaries of a computational domain with a perfectly matched layer (PML) of finite thickness in order to prevent artificially reflected waves from contaminating a numerical simulation. Unfortunately, the PML does not give us an indication about appropriate boundary conditions needed to close the edges of the PML, or how those boundary conditions should be enforced in a numerical setting. Terminating the PML with an inappropriate boundary condition or an unstable numerical boundary procedure can lead to exponential growth in the PML which will eventually destroy the accuracy of a numerical simulation everywhere. In this paper, we analyze the stability and the well-posedness of boundary conditions terminating the PML for the elastic wave equation in first order form. First, we consider a vertical modal PML truncating a two space dimensional computational domain in the horizontal direction. We freeze all coefficients and consider a left half-plane problem with linear boundary conditions terminating the PML. The normal mode analysis is used to study the stability and well-posedness of the resulting initial boundary value problem (IBVP). The result is that any linear well-posed boundary condition yielding an energy estimate for the elastic wave equation, without the PML, will also lead to a well-posed IBVP for the PML. Second, we extend the analysis to the PML corner region where both a horizontal and vertical PML are simultaneously active. The challenge lies in constructing accurate and stable numerical approximations for the PML and the boundary conditions. Third, we develop a high order accurate finite difference approximation of the PML subject to the boundary conditions. To enable accurate and stable numerical boundary treatments for the PML we construct continuous energy estimates in the Laplace space for a one space dimensional problem and two space dimensional PML corner problem. We use summation-by-parts finite difference operators to approximate the spatial derivatives and impose boundary conditions weakly using penalties. In order to ensure numerical stability of the discrete PML, it is necessary to extend the numerical boundary procedure to the auxiliary differential equations. This is crucial for deriving discrete energy estimates analogous to the continuous energy estimates. Numerical experiments are presented corroborating the theoretical results. Moreover, in order to ensure longtime numerical stability, the boundary condition closing the PML, or its corresponding discrete implementation, must be dissipative. Furthermore, the numerical experiments demonstrate the stable and robust treatment of PML corners.

Stable and high order accurate difference methods for the elastic wave equation in discontinuous media

In this paper, we develop a stable and systematic procedure for numerical treatment of elastic waves in discontinuous and layered media. We consider both planar and curved interfaces where media parameters are allowed to be discontinuous. The key feature is the highly accurate and provably stable treatment of interfaces where media discontinuities arise. We discretize in space using high order accurate finite difference schemes that satisfy the summation by parts rule. Conditions at layer interfaces are imposed weakly using penalties. By deriving lower bounds of the penalty strength and constructing discrete energy estimates we prove time stability. We present numerical experiments in two space dimensions to illustrate the usefulness of the proposed method for simulations involving typical interface phenomena in elastic materials. The numerical experiments verify high order accuracy and time stability.

Diagonal-norm summation by parts operators for finite difference approximations of third and fourth derivatives

High-order accurate finite difference operators for third and fourth derivatives are derived. The closures are based on the summation-by-parts (SBP) framework, thereby guaranteeing linear stability. Stability for nonlinear equations that support a convex extension can be achieved if the SBP operators are based on a diagonal norm. The boundary conditions are imposed using a penalty technique. The accuracy and stability properties of the newly derived SBP operators are demonstrated for both linear and nonlinear dispersive wave propagation problems.

Review of summation-by-parts schemes for initial–boundary-value problems

High-order finite difference methods are efficient, easy to program, scale well in multiple dimensions and can be modified locally for various reasons (such as shock treatment for example). The main drawback has been the complicated and sometimes even mysterious stability treatment at boundaries and interfaces required for a stable scheme. The research on summation-by-parts operators and weak boundary conditions during the last 20 years has removed this drawback and now reached a mature state. It is now possible to construct stable and high order accurate multi-block finite difference schemes in a systematic building-block-like manner. In this paper we will review this development, point out the main contributions and speculate about the next lines of research in this area.

High-fidelity numerical solution of the time-dependent Dirac equation

A stable high-order accurate finite difference method for the time-dependent Dirac equation is derived. Grid-convergence studies in 1-D and 3-D corroborate the analysis. The method is applied to time-resolved quantum tunneling where a comparison with the solution to the time-dependent Schrödinger equation in 1-D illustrates the differences between the two equations. In contrast to the conventional tunneling probability decay predicted by the Schrödinger equation, the Dirac equation exhibits Klein tunneling. Solving the time-dependent Dirac equation with a step potential in 3-D reveals that particle spin affects the tunneling process. The observed spin-dependent reflection allows for a new type of spin-selective measurements.

Temporal and spatial evolution characteristics of disturbance wave in a hypersonic boundary layer due to single-frequency entropy disturbance.

By using a high-order accurate finite difference scheme, direct numerical simulation of hypersonic flow over an 8° half-wedge-angle blunt wedge under freestream single-frequency entropy disturbance is conducted; the generation and the temporal and spatial nonlinear evolution of boundary layer disturbance waves are investigated. Results show that, under the freestream single-frequency entropy disturbance, the entropy state of boundary layer is changed sharply and the disturbance waves within a certain frequency range are induced in the boundary layer. Furthermore, the amplitudes of disturbance waves in the period phase are larger than that in the response phase and ablation phase and the frequency range in the boundary layer in the period phase is narrower than that in these two phases. In addition, the mode competition, dominant mode transformation, and disturbance energy transfer exist among different modes both in temporal and in spatial evolution. The mode competition changes the characteristics of nonlinear evolution of the unstable waves in the boundary layer. The development of the most unstable mode along streamwise relies more on the motivation of disturbance waves in the upstream than that of other modes on this motivation.

A perfectly matched layer for the time-dependent wave equation in heterogeneous and layered media

A mathematical analysis of the perfectly matched layer (PML) for the time-dependent wave equation in heterogeneous and layered media is presented. We prove the stability of the PML for discontinuous media with piecewise constant coefficients, and derive energy estimates for discontinuous media with piecewise smooth coefficients. We consider a computational setup consisting of smaller structured subdomains that are discretized using high order accurate finite difference operators for approximating spatial derivatives. The subdomains are then patched together into a global domain by a weak enforcement of interface conditions using penalties. In order to ensure the stability of the discrete PML, it is necessary to transform the interface conditions to include the auxiliary variables. In the discrete setting, the transformed interface conditions are crucial in deriving discrete energy estimates analogous to the continuous energy estimates, thus proving stability and convergence of the numerical method. Finally, we present numerical experiments demonstrating the stability of the PML in a layered medium and high order accuracy of the proposed interface conditions.

Discussion on Direct Numerical Simulation of Large-Scale Structures in High-Reynolds Number Channel Flow by Means of High-Order Accurate Finite Difference Method for Wall-Parallel Directions

Discussion on Direct Numerical Simulation of Large-Scale Structures in High-Reynolds Number Channel Flow by Means of High-Order Accurate Finite Difference Method for Wall-Parallel Directions

Simulation of Dynamic Earthquake Ruptures in Complex Geometries Using High-Order Finite Difference Methods

We develop a stable and high-order accurate finite difference method for problems in earthquake rupture dynamics in complex geometries with multiple faults. The bulk material is an isotropic elastic solid cut by pre-existing fault interfaces that accommodate relative motion of the material on the two sides. The fields across the interfaces are related through friction laws which depend on the sliding velocity, tractions acting on the interface, and state variables which evolve according to ordinary differential equations involving local fields. The method is based on summation-by-parts finite difference operators with irregular geometries handled through coordinate transforms and multi-block meshes. Boundary conditions as well as block interface conditions (whether frictional or otherwise) are enforced weakly through the simultaneous approximation term method, resulting in a provably stable discretization. The theoretical accuracy and stability results are confirmed with the method of manufactured solutions. The practical benefits of the new methodology are illustrated in a simulation of a subduction zone megathrust earthquake, a challenging application problem involving complex free-surface topography, nonplanar faults, and varying material properties.

Convergence Analysis of Weighted Difference Approximations on Piecewise Uniform Grids to a Class of Singularly Perturbed Functional Differential Equations

We consider the numerical approximation of a singularly perturbed time delayed convection diffusion problem on a rectangular domain. Assuming that the coefficients of the differential equation be smooth, we construct and analyze a higher order accurate finite difference method that converges uniformly with respect to the singular perturbation parameter. The method presented is a combination of the central difference spatial discretization on a Shishkin mesh and a weighted difference time discretization on a uniform mesh. A priori explicit bounds on the solution of the problem are established. These bounds on the solution and its derivatives are obtained using a suitable decomposition of the solution into regular and layer components. It is shown that the proposed method is \(L_{2}^{h}\)-stable. The analysis done permits its extension to the case of adaptive meshes which may be used to improve the solution. Numerical examples are presented to demonstrate the effectiveness of the method. The convergence obtained in practical satisfies the theoretical predictions.

Superconvergent Functional Estimates from Summation-By-Parts Finite-Difference Discretizations

Diagonal-norm summation-by-parts (SBP) operators can be used to construct time-stable high-order accurate finite-difference schemes. However, to achieve both stability and accuracy, these operators must use s-order accurate boundary closures when the interior scheme is $2s$-order accurate. The boundary closure limits the solution to $(s+1)$-order global accuracy. Despite this bound on solution accuracy, we show that functional estimates can be constructed that are $2s$-order accurate. This superconvergence requires dual-consistency, which depends on the SBP operators, the boundary condition implementation, and the discretized functional. The theory is developed for scalar hyperbolic and elliptic partial differential equations in one dimension. In higher dimensions, we show that superconvergent functional estimates remain viable in the presence of curvilinear multiblock grids with interfaces. The generality of the theoretical results is demonstrated using a two-dimensional Poisson problem and a nonlinear hyperbolic system—the Euler equations of fluid mechanics.

Moment preserving schemes for Euler equations

A high order accurate finite difference scheme is proposed for one-dimensional Euler equations. In the scheme a set of first three moments of each signal are preserved during the updating. The scheme is one of 5th order in space and 4th order in time. This feature is different from that in typical existing methods in which the use of the first three polynomials results in only 3rd order accuracy in space. The scheme has different features from the existing high order schemes, and the most noticeable are the simultaneous discretization both in space and time, and the use of moments of Riemann invariants instead of primitive physical variables. Numerical examples are given to show the accuracy of the scheme and its robustness for the flows involving shocks.