High-order Finite Difference Method Research Articles

A conservative splitting high-order finite difference method for coupled Gross–Pitaevskii equations in 2D

A conservative splitting high-order finite difference method for coupled Gross–Pitaevskii equations in 2D

Fitted computational method for singularly perturbed convection-diffusion equation with time delay

A uniformly convergent numerical scheme is proposed to solve a singularly perturbed convection-diffusion problem with a large time delay. The diffusion term of the problem is multiplied by a perturbation parameter, ε. For a small ε, the problem exhibits a boundary layer, which makes it challenging to solve it analytically or using standard numerical methods. As a result, the backward Euler scheme is applied in the temporal direction. Non-symmetric finite difference schemes are applied for approximating the first-order derivative terms, and a higher-order finite difference method is applied for approximating the second-order derivative term. Furthermore, an exponential fitting factor is computed and induced in the difference scheme to handle the effect of the small parameter. Using the discrete maximum principle, the stability of the scheme is examined and analyzed. The developed scheme is parameter-uniform with a linear order of convergence in both space and time. To examine the accuracy of the method, two model examples are considered. Further, the boundary layer behavior of the solutions is given graphically.

A robust higher-order finite difference technique for a time-fractional singularly perturbed problem

A higher-order finite difference method is developed to solve the variable coefficients convection–diffusion singularly perturbed problems (SPPs) involving fractional-order time derivative with the order α∈(0,1). The solution to this problem class possesses a typical weak singularity at the initial time t=0 and an exponential boundary layer at the right lateral surface as the perturbation parameter ɛ→0. Alikhanov’s L2−1σ approximation is applied in the temporal direction on a suitable graded mesh, and the spatial variable is discretized on a piecewise uniform Shishkin mesh using a combination of midpoint upwind and central finite difference operators. Stability estimates and the convergence analysis of the fully discrete scheme are provided. It is shown that the fully discrete scheme is uniformly convergent with a rate of O(M−p+N−2(logN)2), where p=min{2,rα}, r is the graded mesh parameter and M,N are the number of mesh points in the time and space direction, respectively. Two numerical examples are taken in counter to confirm the sharpness of the theoretical estimates.

Real-space solution to the electronic structure problem for nearly a million electrons.

We report a Kohn-Sham density functional theory calculation of a system with more than 200 000 atoms and 800 000 electrons using a real-space high-order finite-difference method to investigate the electronic structure of large spherical silicon nanoclusters. Our system of choice was a 20nm large spherical nanocluster with 202 617 silicon atoms and 13 836 hydrogen atoms used to passivate the dangling surface bonds. To speed up the convergence of the eigenspace, we utilized Chebyshev-filtered subspace iteration, and for sparse matrix-vector multiplications, we used blockwise Hilbert space-filling curves, implemented in the PARSEC code. For this calculation, we also replaced our orthonormalization + Rayleigh-Ritz step with a generalized eigenvalue problem step. We utilized all of the 8192 nodes (458 752 processors) on the Frontera machine at the Texas Advanced Computing Center. We achieved two Chebyshev-filtered subspace iterations, yielding a good approximation of the electronic density of states. Our work pushes the limits on the capabilities of the current electronic structure solvers to nearly 106 electrons and demonstrates the potential of the real-space approach to efficiently parallelize large calculations on modern high-performance computing platforms.

A numerical method for two-phase flow with its application to the Kelvin–Helmholtz instability problem

This article uses the artificial compressibility-based high-order finite difference method to simulate the two-phase Kelvin–Helmholtz (KH) instability problem. The equations are based on the mass-conserving Allen–Cahn equation coupled with the incompressible Navier–Stokes equations. One of the advantages of the artificial compressibility approach is that many high-order numerical schemes based on hyperbolic conservation law can be applied. A fifth-order weighted essentially non-oscillatory (WENO) scheme is used for discretizing the convective terms while a standard central finite difference scheme is used for the viscous and surface tension terms. The system of equations is then solved using the Beam-Warming approximate factorization technique. For validation, the effects of both single- and double-mode sinusoidal perturbations on the Kelvin–Helmholtz instability dynamics are analyzed. When there is a single-mode sinusoidal perturbation, the interface roll-up at the center of the domain. Additionally, the role of the surface tension parameter in the instability’s dynamics is investigated. The development of Kelvin–Helmholtz instability is shown to be sensitive to the surface tension value. The analysis of grid convergence is also performed to capture the interface dynamics at varying resolutions. The comparison of the computed results is in good agreement with those from the literature. It is observed that the proposed technique effectively resolves the dynamics of the chaotically distorted interfaces of the Kelvin–Helmholtz instability in two-phase flow.

Solving Lane–Emden equations with boundary conditions of various types using high-order compact finite differences

In this study, a high-order compact finite difference method is used to solve Lane–Emden equations with various boundary conditions. The norm is to use a first-order finite difference scheme to approximate Neumann and Robin boundary conditions, but that compromises the accuracy of the entire scheme. As a result, new higher-order finite difference schemes for approximating Robin boundary conditions are developed in this work. We test the applicability and performance of the method using different examples of Lane–Emden equations. Convergence analysis is provided, and it is consistent with the numerical results. The results are compared with the exact solutions and published results from other methods. The method produces highly accurate results, which are displayed in tables and graphs.

Local-basis difference potentials method for elliptic PDEs in complex geometry

We develop efficient and high-order accurate finite difference methods for elliptic partial differential equations in complex geometry in the difference potentials framework. The main novelty of the developed schemes is the use of local basis functions defined at near-boundary grid points. The use of local basis functions allows unified numerical treatment of (i) explicitly and implicitly defined geometry; (ii) geometry of more complicated shapes, such as those with corners, multi-connected domain, etc; and (iii) different types of boundary conditions. This geometrically flexible approach is complementary to the classical difference potentials method using global basis functions, especially in the case where a large number of global basis functions are needed to resolve the boundary, or where the optimal global basis functions are difficult to obtain. Fast Poisson solvers based on FFT are employed for standard, centered finite difference stencils regardless of the designed order of accuracy. Proofs of convergence of difference potentials in maximum norm are outlined both theoretically and numerically.

A high-order finite difference method with immersed-boundary treatment for fully-nonlinear wave–structure interaction

In order to predict nonlinear wave loading on marine structures, a fully nonlinear higher-order finite difference based potential flow solver with all boundary conditions treated by an immersed boundary method has been developed in this paper. The solver adopts high-order finite difference schemes for the spatial derivatives and a 4th order Runge–Kutta method for time stepping. Test cases of forced oscillation of a cylinder in an infinite fluid domain are first studied, which reveal the advantage of the acceleration potential method in terms of wave load computation. Then a wave generation problem using a piston type wave maker is tested. Special attention is paid to the intersection point between the free surface and the body surface, and a scheme which best meets the accuracy and stability requirements is suggested from several proposals. A novel hyperviscosity filter, which works for both uniform and non-uniform grids, is introduced to stabilize the time-domain solution of the wave maker problem. Finally, a forced heaving cylinder on the free surface is considered, and the nonlinear wave loads on the cylinder are analyzed and compared to benchmark results.

α-, β-phenomena for the supersonic flow around an impulsively started circular cylinder

The vortex dynamics of a supersonic flow past an impulsively started circular cylinder are studied numerically. A well-developed high-order finite difference method is adopted to solve the two-dimensional Navier–Stokes equations. In particular, Ma = 1.1 and Re = 2.2 × 104 are selected to report the vortex evolution process. The results confirm that the α-, β-, sub-α-, and sub-β-phenomena observed in incompressible flow in the separated region do exist in the supersonic regime. In addition, the influence of Mach number and Reynolds number on the vortex is also studied, and the vortex changes are compared at Ma = 1.1, 1.7 and Re = 1 × 104, 2 × 104, 2.5 × 104, 3 × 104, 5 × 104, and 1 × 105.

On construction of shock-capturing boundary closures for high-order finite difference method

High-order finite difference method is widely used in the field of computational fluid dynamics. Since it can be implemented in a dimension-by-dimension fashion, it is much more efficient than some other methods in higher dimensions. However, although interior shock-capturing schemes are mature, shock-capturing boundary closures are only in their infancy due to the issue of insufficient stencils near boundary. To address this issue, we focus on a fifth-order multi-resolution weighted compact nonlinear scheme. Firstly, to enhance the shock-capturing capability of the interior scheme, we propose a new technique by adding a free parameter to the formulae of smoothness indicators. Then we follow the idea of the nonlinear spectral difference method to construct shock-capturing boundary closures when the interior scheme cannot be used. To validate the proposed method, several benchmark numerical examples are presented. The results show that the shock-capturing capability is effectively improved and the proposed nonlinear boundary closures work well for the cases with shock transmission or shock near boundary.

Simulations of interaction between flow and a thin elastic beam using CIP-based model with adaptive mesh

The simulation of interaction between flow and the elastic structure is a complicated fluid-solid interaction problem, involving the structure's large deformation and fluid interface variation. In this paper, an adaptive CIP (constrained interpolation profile)-based numerical model has been presented to simulate this problem. In this CIP-based model, the constrained interpolation profile method, a high-order finite difference method, is employed as the base flow solver and THINC/SW (Tangent of Hyperbola for Interface Capturing with Slope Weighting) method, a volume of fluid method, is applied to capture the flow interface. IBM (immersed boundary method) method is used to treat the fluid-solid interface. The body deformation is calculated by coupled FDM-FEM (finite difference method and finite element method) method. Adaptive mesh is generated by Afivo (Adaptive Finite Volume quadtree/Octree mesh) framework. The cases of flow past a square and water entry of a symmetric solid wedge are carried out to validate the reliability of the flow solver in this present model. Cantilever beam deformation case is used to validate the accuracy of the solid solver in the present model. Then cases of a thin elastic beam in sloshing tank and dam-break flow are carried out by using this model. After that, the case of regular wave past a membrane breakwater is studied. Our numerical results are compared with previous experimental and numerical studies, and good agreements are obtained. It is demonstrated that this present numerical model has good applicability to simulate the interaction between flow and a thin elastic beam.

A Compact High-Order Finite-Difference Method with Optimized Coefficients for 2D Acoustic Wave Equation

High-precision finite difference (FD) wavefield simulation is one of the key steps for the successful implementation of full-waveform inversion and reverse time migration. Most explicit FD schemes for solving seismic wave equations are not compact, which leads to difficulty and low efficiency in boundary condition treatment. Firstly, we review a family of tridiagonal compact FD (CFD) schemes of various orders and derive the corresponding optimization schemes by minimizing the error between the true and numerical wavenumber. Then, the optimized CFD (OCFD) schemes and a second-order central FD scheme are used to approximate the spatial and temporal derivatives of the 2D acoustic wave equation, respectively. The accuracy curves display that the CFD schemes are superior to the central FD schemes of the same order, and the OCFD schemes outperform the CFD schemes in certain wavenumber ranges. The dispersion analysis and a homogeneous model test indicate that increasing the upper limit of the integral function helps to reduce the spatial error but is not conducive to ensuring temporal accuracy. Furthermore, we examine the accuracy of the OCFD schemes in the wavefield modeling of complex structures using a Marmousi model. The results demonstrate that the OCFD4 schemes are capable of providing a more accurate wavefield than the CFD4 scheme when the upper limit of the integral function is 0.5π and 0.75π.

Wave field characteristics and transmission processing of seismic waves in coal seams with non-developing fundamental mode channel waves

Fine detection of geological structures in the working face of mines can provide a strong geological guarantee for safe and efficient production. The channel wave seismic method is an effective geological guarantee technology; however, it is only effective for coal seams wherein channel waves develop well. It does not yield acceptable results for coal seams in which channel waves develop poorly or fail to develop. To address this problem, the Kelvin–Voigt viscoelastic model and the staggered grid high-order finite difference method were applied to study the propagation and attenuation characteristics of the seismic wave field in viscoelastic coal seams. Moreover, the causes for the non-development of channel waves were explored. Numerical simulation and a field experiment on transmitted channel waves were performed for fault structure detection in coal seams. The effects of refracted longitudinal waves, multiple shear horizontal (SH) waves, and first-order Rayleigh channel waves on fault structure detection in models with different fault throws were investigated. The results indicate that the coal seam transverse wave quality factor Qs has a significant influence on the absorption attenuation of Love channel waves and fundamental mode Rayleigh channel waves. However, it has minimal influence on refracted longitudinal waves, first-order Rayleigh channel waves, and multiple SH waves. When the fault throw is greater than or equal to the coal seam thickness, the first-order Rayleigh channel waves and multiple SH waves can accurately reflect the fault location and extension shape, and the refracted P waves can accurately locate the faults but cannot accurately reflect the fault extension shape. When the fault throw is half the coal seam thickness, the refracted P waves and multiple SH waves have a poor detection effect, but the first-order Rayleigh channel waves can accurately reflect the fault location and extension shape. For working faces in which fundamental mode channel waves develop poorly, first-order Rayleigh channel waves can replace the fundamental mode channel waves in seismic methods to locate geological structures in coal seam working faces.

Efficient numerical schemes for multidimensional population balance models

Multidimensional population balance models (PBMs) describe chemical and biological processes having a distribution over two or more intrinsic properties (such as size and age, or two independent spatial variables). The incorporation of additional intrinsic variables into a PBM improves its descriptive capability and can be necessary to capture specific features of interest. As most PBMs of interest cannot be solved analytically, computationally expensive high-order finite difference or finite volume methods are frequently used to obtain an accurate numerical solution. We propose a finite difference scheme based on operator splitting and solving each sub-problem at the limit of numerical stability that achieves a discretization error that is zero for certain classes of PBMs and low enough to be acceptable for other classes. In conjunction to employing specially constructed meshes and variable transformations, the scheme exploits the commutative property of the differential operators present in many classes of PBMs. The scheme has very low computational cost — potentially as low as just memory reallocation. Multiple case studies demonstrate the performance of the proposed scheme.

Effect of anisotropic resistance characteristic on boundary-layer transitional flow

It is observed that the feather surface exhibits anisotropic resistances for the streamwise and spanwise flows. To obtain a qualitative understanding about the effect of this anisotropic resistance feature of surface on the boundary-layer transitional flow over a flat plate, a simple phenomenological model for the anisotropic resistance is established in this paper. By means of the large eddy simulation (LES) with high-order accurate finite difference method, the numerical investigations are conducted. The numerical results show that with the spanwise resistance hindering the formation of vortexes, the transition from laminar flow to turbulent flow can be delayed, and turbulence is weakened when the flow becomes fully turbulent, which leads to significant drag reduction for the plate. On the contrary, the streamwise resistance renders the flow less stable, which leads to the earlier transition and enhances turbulence in the turbulent region, causing a drag increase for the plate. Thus, it is indicated that a surface with large resistance for spanwise flow and small resistance for streamwise flow can achieve significant drag reduction. The present results highlight the anisotropic resistance characteristic near the feather surface for drag reduction, and shed a light on the study of bird’s efficient flight.

High-Order Compact Finite Difference Methods for Solving the High-Dimensional Helmholtz Equations

Abstract In this paper, the sixth-order compact finite difference schemes for solving two-dimensional (2D) and three-dimensional (3D) Helmholtz equations are proposed. Firstly, the sixth-order compact difference operators for the second-order derivatives are applied to approximate the Laplace operator. Meanwhile, with the original differential equation, the sixth-order compact difference schemes are proposed. However, the truncation errors of the proposed scheme obviously depend on the unknowns, source function and wavenumber. Thus, we correct the truncation error of the sixth-order compact scheme to obtain an improved sixth-order compact scheme that is more accurate. Theoretically, the convergence and stability of the present improved method are proved. Finally, numerical tests verify that the improved schemes are more accurate.

Highly Accurate Compact Finite Difference Schemes for Two-Point Boundary Value Problems with Robin Boundary Conditions

In this study, a high-order compact finite difference method is used to solve boundary value problems with Robin boundary conditions. The norm is to use a first-order finite difference scheme to approximate Neumann and Robin boundary conditions, but that compromises the accuracy of the entire scheme. As a result, new higher-order finite difference schemes for approximating Robin boundary conditions are developed in this work. Six examples for testing the applicability and performance of the method are considered. Convergence analysis is provided, and it is consistent with the numerical results. The results are compared with the exact solutions and published results from other methods. The method produces highly accurate results, which are displayed in tables and graphs.

Convergence Rate of the High-Order Finite Difference Method for Option Pricing in a Markov Regime-Switching Jump-Diffusion Model

The high-order finite difference method for option pricing is one of the most popular numerical algorithms. Therefore, it is of great significance to study its convergence rate. Based on the relationship between this algorithm and the trinomial tree method, as well as the definition of local remainder estimation, a strict mathematical proof is derived for the convergence rate of the high-order finite difference method for option pricing in a Markov regime-switching jump-diffusion model. The theoretical result shows that the convergence rate of this algorithm is O(Δτ) . Moreover, the results also hold in the case of Brownian motion and jump-diffusion models that are specialized forms of the given model.

A High Order Finite Difference Method for the Elastic Wave Equation in Bounded Domains with Nonconforming Interfaces

A High Order Finite Difference Method for the Elastic Wave Equation in Bounded Domains with Nonconforming Interfaces

A higher-order finite difference method for two-dimensional singularly perturbed reaction-diffusion with source-term-discontinuous problem

This paper considers a two-dimensional singularly perturbed reaction-diffusion equation with a discontinuous source term. Due to this discontinuity, interior, corner, and boundary layers appear in the solution for adequately small values of the perturbation parameter ϵ. To achieve a decent estimate of the solution, we construct a numerical approach adopting an efficient hybrid finite difference method that includes a proper layer adapted piece-wise uniform Shishkin mesh. Further, we prove that the hybrid finite difference method is almost second-order uniformly convergent with respect to the perturbation parameter. We have implemented our method to test examples. Numerical results are verifying the theoretical results.