High-order Finite Difference Method Research Articles

Exponential high-order compact finite difference method for convection-dominated diffusion problems on nonuniform grids

In this article, an exponential high-order compact (EHOC) difference method is introduced for solving the convection-dominated convection diffusion problems on nonuniform grid without any coordinate transformation from the physical space to the computational space. The derived EHOC schemes on nonuniform grid can not only preserve nonoscillation property and yield high accuracy approximation solutions, but also efficiently handle the convection diffusion problems with boundary layers by employing a flexible discretized grid that can be adapted to the singularity in the domain. To demonstrate the performances in computational accuracy, efficiency, and stability of the proposed EHOC schemes, several problems with boundary or interior layers where sharp gradients may appear are solved numerically. The results obtained by the present EHOC schemes on nonuniform grid are compared with analytical solutions and those using the EHOC on uniform grid and other methods. The present EHOC schemes on nonuniform grid provide excellent results for all test problems. It is shown that the proposed EHOC schemes on nonuniform grid are accurate, efficient, and stable, as well as have the advantage of better-scale resolution for convection-dominated diffusion problems.

Accurate gradient preserved method for solving heat conduction equations in double layers

Analyzing heat transfer in layered structures is important for the design and operation of devices and the optimization of thermal processing of materials. Existing numerical methods dealing with layered structures if using only three grid points across the interface usually provide only a second-order truncation error. Obtaining a numerical scheme using three grid points across the interface so that the overall numerical scheme is unconditionally stable and convergent with higher-order accuracy is mathematically challenging. In this study, we develop a higher-order accurate finite difference method using three grid points across the interface by preserving the first-order derivative, ux, in the interfacial condition and/or the boundary condition. As such, when the compact fourth-order accurate Padé scheme is used at the interior points, the overall scheme maintains to be higher-order accurate. Stability and convergence of the scheme are analyzed. Finally, four examples are given to test the obtained numerical method. Results show that the convergence order in space is close to 4.0, which coincides with the theoretical analysis.

Numerical simulation of the seismic wave propagation and fluid pressure in complex porous media at the mesoscopic scale

In this study, we apply an improved high-order staggered-grid finite difference method to model seismic wave propagation in complex heterogeneous porous media with cracks. Upon changing the mesh step and applying a broadband wavelet, we obtain numerical simulation results at the mesoscopic scale. As a fundamental result of applying Biot’s theory, three types of waves are also found. We discuss the relationship between the shear wave splitting parameter and the crack density of an equivalent medium in detail. Based on Darcy’s law, the finite element method is used to investigate the distribution of the fluid pressure field in porous media at the mesoscopic scale. A flat ellipse is used as an approximation of a crack. According to the simulation results, we further study the influence of variations in the crack angle on the distribution of the fluid pressure. Moreover, the formulas derived from the Fourier transform are used to calculate the inverse attenuation quality factors and phase velocities of the fast compressional wave, slow compressional wave and shear wave. The effects of the porosity and elastic modulus of an equivalent medium on wave attenuation are also taken into account in this study.

Numerical study of air compressibility effects in breaking wave impacts using a CIP-based model

Abatract The prediction of breaking wave loads is challenging in ocean engineering due to air entrapment and compressibility. In order to investigate air compressibility effects in different modes of breaking wave impact, a two-phase compressible flow model is proposed. A high-order finite difference method, Constrained Interpolation Profile (CIP) method is employed to discretize the Navier-Stokes equations, and the interface is captured by Tangent of Hyperbola for Interface Capturing (THINC) scheme. The present compressible model is validated against several tests: 1-D gravity-induced liquid piston motion, idealized free drop of water patch impact and dam-break experiments. Special attention of this study is paid to breaking wave impacts on a vertical wall in a liquid sloshing tank. For air entrapped impact, the results of CIP-based compressible model are in good agreement with experimental data. It is shown that the air compressibility leads to high-frequency oscillation of the impact pressure, which cannot be reproduced by the incompressible model. It is found that the oscillation of the impact pressure is synchronized with the compression and expansion of the entrapped air pocket, and the oscillation frequency is related to the size of air pocket. For air-water mixing impact, the oscillation of the pressure may be more intense.

Interaction of wake disturbance with compressible transitional boundary layers in a low-pressure turbine cascade under rotor-stator interaction

Transition-process-resolving numerical simulations of compressible transitional flows subjected to rotor-stator interaction in a low-pressure turbine cascade were performed. Cases of small and large Strouhal numbers (St) of blade passage, and with and without isotropic free-stream turbulence (IFST) are compared. Cascade geometry is T106, and the computations are conducted at its design point condition. A higher-order finite difference method based on sixth-order compact difference and tenth-order filtering is employed. Variation in St and the intensity of IFST changes the frequency of turbulence inflow, the convection and propagation of turbulence, and the interaction of disturbance with boundary layers in the blade passage. The changes result in different pressure fluctuation in the passage. Because pressure fluctuation can convey information away from the passage, it is expected to be used for diagnosing boundary layer states around the blade indirectly away from the passage.

A variable high-order shock-capturing finite difference method with GP-WENO

We present a new finite difference shock-capturing scheme for hyperbolic equations on static uniform grids. The method provides selectable high-order accuracy by employing a kernel-based Gaussian Process (GP) data prediction method which is an extension of the GP high-order method originally introduced in a finite volume framework by the same authors. The method interpolates Riemann states to high order, replacing the conventional polynomial interpolations with polynomial-free GP-based interpolations. For shocks and discontinuities, this GP interpolation scheme uses a nonlinear shock handling strategy similar to Weighted Essentially Non-Oscillatory (WENO), with a novelty consisting in the fact that nonlinear smoothness indicators are formulated in terms of the Gaussian likelihood of the local stencil data, replacing the conventional L2-type smoothness indicators of the original WENO method. We demonstrate that these GP-based smoothness indicators play a key role in the new algorithm, providing significant improvements in delivering high – and selectable – order accuracy in smooth flows, while successfully delivering non-oscillatory solution behavior in discontinuous flows.

Land use Detection in Nusajaya using Higher-Order Modified Geodesic Active Contour Model

Urban development is a global phenomenon. In Johor, especially Nusajaya is one of the most rapidly developing cities. This is due to the increase of land demand and population growth. Moreover, land-use changes are considered to be one of the major components of current environmental monitoring strategies. In this context, image segmentation and mathematical model offers essential tools that can be used to analyze land use detection. The image segmentation process is known as the most important and difficult task in image analysis. Nonlinear fourth-order models had shown to have a good achievement in recovering smooth regions. Therefore, these motivate us to propose a fourth-order modified geodesic active contour (GAC) model. In the proposed model, a modified signed pressure force (SPF) function has been defined to segment the inhomogeneous satellite images. The simulations of the fourth-order modified GAC model through some numerical methods based on the higher-order finite difference method (FDM) have been illustrated. Matlab R2015a software in Windows 7 Ultimate on Intel (R) Core (TM) i5-3230M @ 2.60GHz CPU with 8 GB RAM has been considered as a computational platform for the simulation. Qualitative and quantitative differences between modified SPF functions and other SPF functions have been shown as a comparison. Hence land use detection is very useful for local governments and urban planners to enhance the future sustainable development plans of Nusajaya.

High-Order Accurate Conservative Finite Difference Methods for Vlasov Equations in 2D+2V

This manuscript discusses discretization of the Vlasov--Poisson system in 2D+2V phase space using high-order accurate conservative finite difference algorithms. One challenge confronting direct kinetic simulation is the significant computational cost associated with the high-dimensional phase space description. In the present work we advocate the use of high-order accurate schemes as a mechanism to reduce the computational cost required to deliver a given level of error in the computed solution. We pursue a discretely conservative finite difference formulation of the governing equations, and discuss fourth- and sixth-order accurate schemes. In addition, we employ a minimally dissipative nonlinear scheme based on the well-known WENO (weighted essentially nonoscillatory) approach. Verification of the full formulation is performed using the method of manufactured solutions. Results are also presented for the physically relevant scenarios of Landau damping, and growth of transverse instabilities from an imposed plane wave.

High-Order Finite Difference Method for Helmholtz Equation in Polar Coordinates

We present a fourth-order finite difference scheme for the Helmholtz equation in polar coordinates. We employ the finite difference format in the interior of the region and derive a nine-point fourth-order scheme. Specially, ghost points outside the region are applied to obtain the approximation for the Neumann boundary condition. We obtain the matrix form of the linear system and the sparsity of the coefficient matrix is favorable for the computation of the Helmholtz equation. The feasibility and accuracy of the method are validated by two test examples which have exact solutions.

Three-dimensional direct numerical simulation of infrasound propagation in the Earth’s atmosphere

A direct numerical simulation of the three-dimensional unsteady compressible Navier–Stokes equations is performed to investigate the infrasonic field generated in a realistic atmosphere by an explosive source placed at ground level. To this end, a high-order finite-difference method originally developed for aeroacoustic applications is employed. The maximum overpressure and the main frequency of the signal recorded at 4 km distance from the source location are about 4000 Pa and 0.2 Hz, respectively. The atmosphere is parametrized as a vertically stratified medium, constructed by specifying vertical profiles of the temperature and the horizontal wind which reproduce measurements. The computation is carried out up to 140 km altitude and 450 km range. The goal of the present paper is twofold. On the one hand, the feasibility of using a direct numerical simulation of the three-dimensional fluid dynamic equations for the detailed description of long-range propagation in the atmosphere is proven. On the other hand, a physical analysis of the infrasonic field is realized. In particular, great attention is directed towards some important phenomena which are not taken into account or not well predicted by classical propagation models. To begin with, the present study clearly demonstrates that the weakly nonlinear ray theory may lead to an incorrect evaluation of the waveform distortion of high-amplitude waves propagating towards the lower thermosphere. In addition, signals recorded in the shadow zones are investigated. In this regard, the influence on the acoustic field of temperature and wind inhomogeneities of length scale comparable with the acoustic wavelength is analysed. The role of diffraction at the thermospheric caustic is finally examined and it is pointed out that the amplitude of the source may have a strong impact on the length of the shadow zone.

New very high-order upwind multi-layer compact (MLC) schemes with spectral-like resolution for flow simulations

Numerical simulations of multi-scale flow problems such as hypersonic boundary layer transition, turbulent flows, computational aeroacoustics and various other flow problems with complex physics require high-order methods with high spectral resolutions. For instance, the receptivity mechanisms in the hypersonic boundary layer are the resonant interactions between forcing waves and boundary-layer waves, and the complex wave interactions are difficult to be accurately predicted by conventional low-order numerical methods. High-order methods, which are robust and accurate in resolving a wide range of time and length scales, are required. Currently, the high-order finite difference methods for simulations of hypersonic flows are usually upwind schemes or compact schemes with fifth-order accuracy or lower [1]. The objective of this paper is to develop and analyze a new very high-order numerical scheme with the spectral-like resolution for flow simulations on structured grids, with focus on smooth flow problems involving multiple scales. Specifically, a new upwind multi-layer compact (MLC) scheme with spectral-like resolution up to seventh order is derived in a finite difference framework. By using the ‘multi-layer’ idea, which introduces first derivatives into the MLC schemes and approximates the second derivatives, the resolution of the MLC schemes can be significantly improved within a compact grid stencil. The auxiliary equations are required and they are the only nontrivial equations, which contributes to good computational efficiency. In addition, the upwind MLC schemes are derived based on the idea of constructing upwind schemes on centered stencils with adjustable parameters to control the dissipation. Fourier analysis is performed to show that the MLC schemes have small dissipation and dispersion in a very wide range of wavenumbers in both one- and two-dimensional cases, and the anisotropic error is much smaller than conventional finite difference methods in the two-dimensional case. Furthermore, the stability analysis with matrix method shows that high-order boundary closure schemes are stable because of compactness of the stencils. The accuracies and rates of convergence of the new schemes are validated by numerical experiments of the linear advection equation, the nonlinear Euler equations, and the Navier–Stokes equations. The numerical results show that good computational efficiency, very high-order accuracies, and high spectral resolutions especially on coarse meshes can be attained with the MLC schemes. Overall, the MLC scheme has the properties of simple formulations, high-order accuracies, spectral-like resolutions, and compact stencils, and it is suitable for accurate simulation of smooth multi-scale flows with complex physics.

High Accuracy Split-Step Finite Difference Method for Schrödinger-KdV Equations**Supported by the National Natural Science Foundation of China under Grant No. 11571181

In this article, two split-step finite difference methods for Schrödinger-KdV equations are formulated and investigated. The main features of our methods are based on: (i) The applications of split-step technique for Schrödinger-like equation in time. (ii) The utilizations of high-order finite difference method for KdV-like equation in spatial discretization. (iii) Our methods are of spectral-like accuracy in space and can be realized by fast Fourier transform efficiently. Numerical experiments are conducted to illustrate the efficiency and accuracy of our numerical methods.

A high-order L2-compact difference method for Caputo-type time-fractional sub-diffusion equations with variable coefficients

A high-order compact finite difference method is proposed for solving a class of time-fractional sub-diffusion equations. The diffusion coefficient of the equation may be spatially variable and the time-fractional derivative is in the Caputo sense with the order α ∈ (0, 1). The Caputo time-fractional derivative is discretized by a (3−α) th-order numerical formula (called the L2 formula here) which is constructed by piecewise quadratic interpolating polynomials but does not require any sub-stepping scheme for the approximation at the first-time level. The variable coefficient spatial differential operator is approximated by a fourth-order compact finite difference operator. By developing a technique of discrete energy analysis, a full theoretical analysis of the stability and convergence of the method is carried out for the general case of variable coefficient and for all α ∈ (0, 1). The optimal error estimate is obtained in the L2 norm and shows that the proposed method has the temporal (3−α) th-order accuracy and the spatial fourth-order accuracy. Further approximations are also considered for enlarging the applicability of the method while preserving its high-order accuracy. Applications are given to three model problems, and numerical results are presented to demonstrate the theoretical analysis results.

Pseudo-impulsive solutions of the forward-speed diffraction problem using a high-order finite-difference method

This paper considers pseudo-impulsive numerical solutions to the forward-speed diffraction problem, as derived from classical linearized potential flow theory. Both head- and following-seas cases are treated. Fourth-order finite-difference approximations are applied on overlapping, boundary-fitted grids to obtain solutions using both the Neumann-Kelvin and the double-body flow linearizations of the problem. A method for computing the pseudo-impulsive incident wave forcing in finite water depth using the Fast Fourier Transform (FFT) is presented. The pseudo-impulsive scattering solution is then Fourier transformed into the frequency domain to obtain the wave excitation forces and the body motion response. The calculations are validated against reference solutions for a submerged circular cylinder and a submerged sphere. Calculations are also made for a modern bulk carrier, showing good agreement with experimental measurements.

Efficient compact finite difference method for variable coefficient fractional sub-diffusion equations with nonhomogeneous Neumann boundary conditions in conservative form

In this paper, a high-order compact finite difference scheme is derived for the time fractional sub-diffusion equation with nonhomogeneous Neumann boundary conditions in conservative form. Based on the L2- $$1_{\sigma }$$ approximation formula of the time fractional derivative, a high-order efficient compact finite difference method is developed. The unconditional stability of the resulting scheme and its convergence of second order in time and fourth order in space are rigourously proved by a discrete energy analysis method. The theoretical results are illustrated by examples.

Effect of a Roughness Element on the Receptivity of a Hypersonic Boundary Layer over a Blunt Cone Due to Pulse Entropy Disturbance with a Single Frequency.

A high-order finite difference method was used to simulate the hypersonic flow field over a blunt cone with different height roughness elements. The unsteady flow field induced by pulse disturbances was analyzed and compared with that under continuous disturbances. The temporal and spatial evolution characteristics of disturbances in the boundary layer were investigated and the propagation of different disturbance modes in the boundary layer was researched through the fast Fourier transform (FFT) method. The effect of the roughness element on the receptivity characteristic of the hypersonic boundary layer under pulse entropy disturbances was explored. The results showed that the different mode disturbances near roughness in the boundary layer were enlarged in the upstream half of the roughness element and suppressed in the downstream half. However, the effect of roughness weakened gradually as the disturbance frequency increased in the boundary layer. A phenomenon of mode competition in the downstream region of the roughness element exited. As the disturbances propagated downstream, the fundamental mode gradually became the dominant mode. A certain promotion effect on the mode competition was induced by the roughness element and the effect was enhanced with the increase in the roughness element height.

Rock moduli estimation of inhomogeneous two-phase media with finite difference modeling algorithm

Digital rock physics (DRP) builds a bridge between pore-scale physical processes and the macroscopic physical properties of rock. Its key paradigm is to image and digitize the pore space and mineral matrix of rock, then to simulate the response of various physical fields and calculate the equivalent elastic parameters of rocks through digital rocks and mathematical methods. In this paper, a new approach is proposed to estimate the rock moduli of two-phase media with the staggered-grid high-order finite difference method (FDM) of the Biot theory based on DRP. This new method not only takes into consideration the impact of fluid on elastic wave propagation in two-phase media, but it is also easy to understand and implement, improving the calculation accuracy, requiring less memory and improving on the weaknesses of conventional rock physics experiments which are time consuming and expensive. In order to estimate the rock moduli, we establish two models, and the digital rock sample is embedded in one of those. Using this method, it is possible to model the dynamic wave propagation and measure the time delay of the peak amplitude caused by the inhomogeneous structure of the digital rock sample, with the receivers set at the bottom of the two models. The time delay allows us to estimate the effective velocity of both compressional and shear waves, and therefore calculate the rock moduli. Additionally, comparison between the numerical simulated results obtained through this method and experimental results indicates that they agree well. Comparison with the numerical simulated results obtained via another method tests and verifies the accuracy and feasibility of the new method. Also, the equivalence conditions between this new method and the various rock physics models are inferred.

An Improved High Order Finite Difference Method for Non-conforming Grid Interfaces for the Wave Equation

This paper presents an extension of a recently developed high order finite difference method for the wave equation on a grid with non-conforming interfaces. The stability proof of the existing methods relies on the interpolation operators being norm-contracting, which is satisfied by the second and fourth order operators, but not by the sixth order operator. We construct new penalty terms to impose interface conditions such that the stability proof does not require the norm-contracting condition. As a consequence, the sixth order accurate scheme is also provably stable. Numerical experiments demonstrate the improved stability and accuracy property.

A high-order numerical method for solving the 2D fourth-order reaction-diffusion equation

In the present work, orthogonal spline collocation (OSC) method with convergence order O(τ3−α + hr+ 1) is proposed for the two-dimensional (2D) fourth-order fractional reaction-diffusion equation, where τ, h, r, and α are the time-step size, space size, polynomial degree of space, and the order of the time-fractional derivative (0 < α < 1), respectively. The method is based on applying a high-order finite difference method (FDM) to approximate the time Caputo fractional derivative and employing OSC method to approximate the spatial fourth-order derivative. Using the argument developed recently by Lv and Xu (SIAM J. Sci. Comput. 38, A2699–A2724, 2016) and mathematical induction method, the optimal error estimates of proposed fully discrete OSC method are proved in detail. Then, the theoretical analysis is validated by a number of numerical experiments. To the best of our knowledge, this is the first proof on the error estimates of high-order numerical method with convergence order O(τ3−α + hr+ 1) for the 2D fourth-order fractional equation.

A robust method of computing finite difference coefficients based on Vandermonde matrix

When the finite difference (FD) method is employed to simulate the wave propagation, high-order FD method is preferred in order to achieve better accuracy. However, if the order of FD scheme is high enough, the coefficient matrix of the formula for calculating finite difference coefficients is close to be singular. In this case, when the FD coefficients are computed by matrix inverse operator of MATLAB, inaccuracy can be produced. In order to overcome this problem, we have suggested an algorithm based on Vandermonde matrix in this paper. After specified mathematical transformation, the coefficient matrix is transformed into a Vandermonde matrix. Then the FD coefficients of high-order FD method can be computed by the algorithm of Vandermonde matrix, which prevents the inverse of the singular matrix. The dispersion analysis and numerical results of a homogeneous elastic model and a geophysical model of oil and gas reservoir demonstrate that the algorithm based on Vandermonde matrix has better accuracy compared with matrix inverse operator of MATLAB.