High-order Interpolation Research Articles

A Split-Step Algorithm for Effectively Suppressing the Numerical Dispersion for 3D Seismic Propagation Modeling

In this article, we propose a 3D split-step algorithm (SSA) for seismic- wave simulation. We first transform the wave equations in 3D anisotropic media into a system of first-order partial differential equations with respect to time t. Then we use the multidimensional high-order interpolation method to approximate the high-order spatial derivatives, so that we obtain a system of semidiscrete ordinary differential equations (ODEs). Finally, the third-order implicit Adams method and truncated dif- ferentiator series method are applied to solve the semidiscrete ODEs. We provide the theoretical study on the properties of the 3D SSA, such as stability criteria, theoretical error, numerical error, numerical dispersion, and computational efficiency. We also compare some seismic modeling results of this method against those of some high-order finite-difference schemes. Theoretical analysis and numerical tests show that the 3D SSA is third-order accurate in time and fourth-order accurate in space. However, its computational costs and memory requirements are much less than those of the fourth-order Lax-Wendroff correction method and the fourth-order staggered- grid method. Using a multilayer elastic model with large velocity contrasts and free surface, we compare the result of the 3D SSA with that of the discrete-wavenumber method. We also present the synthetic seismograms in the 3D three-layer isotropic medium, the wave-field snapshots in the 3D two-layer medium, and the 3D trans- versely isotropic medium with a vertical symmetry axis. All these promising numer- ical results illustrate that the 3D SSA can suppress effectively the numerical dispersion caused by discretizing the wave equations when too few sampling points per mini- mum wavelength are used or when models have large velocity contrasts between adjacent layers, further resulting in both increasing the computational efficiency and saving the storage space when big grids are used.

Comparison of improved meshless interpolation schemes for SPH method and accuracy analysis

In the smoothed particle hydrodynamics (SPH) method, a meshless interpolation scheme is needed for the unknown function in order to discretize the governing equation. A particle approximation method has so far been used for this purpose. Traditional particle interpolation (TPI) is simple and easy to do, but its low accuracy has become an obstacle to its wider application. This can be seen in the cases of particle disorder arrangements and derivative calculations. There are many different methods to improve accuracy, with the moving least square (MLS) method one of the most important meshless interpolation methods. Unfortunately, it requires complex matrix computing and so is quite time-consuming. The authors developed a simpler scheme, called higher-order particle interpolation (HPI). This scheme can get more accurate derivatives than the MLS method, and its function value and derivatives can be obtained simultaneously. Although this scheme was developed for the SPH method, it has been found useful for other meshless methods.

Explicit Multi-Degree Reduction of Said-Bézier Generalized Ball Curves

给出了计算Said-Bézier型广义Ball曲线(SBGB曲线)在L₂范数下保持端点约束的一种最佳降多阶算法.基于SBGB基函数、幂基函数和Jacobi基函数之间的相互转换关系,得到了SBGB基函数和Jacobi基函数之间的显式转换矩阵;进一步利用Jacobi基的正交性和上述转换矩阵的逆矩阵,导出了SBGB曲线在L₂范数下的显式约束降多阶算法.此算法蕴含了Said-Ball曲线、Bézier曲线以及位置介于这两类曲线之间的一大类参数曲线的相应降多阶算法.证明了这是一种可以预报最佳误差且满足端点高阶约束的一次性降多阶算法.最后用数值实例说明了算法的正确性和优越性.

Compact fourth-order finite volume method for numerical solutions of Navier–Stokes equations on staggered grids

The development of a compact fourth-order finite volume method for solutions of the Navier–Stokes equations on staggered grids is presented. A special attention is given to the conservation laws on momentum control volumes. A higher-order divergence-free interpolation for convective velocities is developed which ensures a perfect conservation of mass and momentum on momentum control volumes. Three forms of the nonlinear correction for staggered grids are proposed and studied. The accuracy of each approximation is assessed comparatively in Fourier space. The importance of higher-order approximations of pressure is discussed and numerically demonstrated. Fourth-order accuracy of the complete scheme is illustrated by the doubly-periodic shear layer and the instability of plane-channel flow. The efficiency of the scheme is demonstrated by a grid dependency study of turbulent channel flows by means of direct numerical simulations. The proposed scheme is highly accurate and efficient. At the same level of accuracy, the fourth-order scheme can be ten times faster than the second-order counterpart. This gain in efficiency can be spent on a higher resolution for more accurate solutions at a lower cost.

A new approach of a limiting process for multi-dimensional flows

An enhanced Multi-dimensional Limiting Process (e-MLP) is developed for the accurate and efficient computation of multi-dimensional flows based on the Multi-dimensional Limiting Process (MLP). The new limiting process includes a distinguishing step and an enhanced multi-dimensional limiting process. First, the distinguishing step, which is independent of high order interpolation and flux evaluation, is newly introduced. It performs a multi-dimensional search of a discontinuity. The entire computational domain is then divided into continuous, linear discontinuous and nonlinear discontinuous regions. Second, limiting functions are appropriately switched according to the type of each region; in a continuous region, there are no limiting processes and only higher order accurate interpolation is performed. In linear discontinuous and nonlinear discontinuous regions, TVD criterion and MLP limiter are respectively used to remove oscillation. Hence, e-MLP has a number of advantages, as it incorporates useful features of MLP limiter such as multi-dimensional monotonicity and straightforward extensionality to higher order interpolation. It is applicable to local extrema and prevents excessive damping in a linear discontinuous region through application of appropriate limiting criteria. It is efficient because a limiting function is applied only to a discontinuous region. In addition, it is robust against shock instability due to the strict distinction of the computational domain and the use of regional information in a flux scheme as well as a high order interpolation scheme. This new limiting process was applied to numerous test cases including one-dimensional shock/sine wave interaction problem, oblique stationary contact discontinuity, isentropic vortex flow, high speed flow in a blunt body, planar shock/density bubble interaction, shock wave/vortex interaction and, particularly, magnetohydrodynamic (MHD) cloud-shock interaction problems. Through these tests, it was verified that e-MLP substantially enhances the accuracy and efficiency with both continuous and discontinuous multi-dimensional flows.

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

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

Spectral element methods on unstructured meshes: which interpolation points?

In the field of spectral element approximations, the interpolation points can be chosen on the basis of different criteria, going from the minimization of the Lebesgue constant to the simplicity of the point generation procedure. In the present paper, we summarize some recent nodal distributions for a high order interpolation in the triangle. We then adopt these points as approximation points for the numerical solution of an elliptic partial differential equation on an unstructured simplicial mesh. The L 2-norm of the approximation error is then analyzed for a model problem.

Simulation of acoustic wavefields in heterogeneous media: A robust method for automatic suppression of numerical dispersion

Numerical dispersion limits the application of numerical simulation methods for solving the acoustic wave equation in large-scale computation. The nearly analytic discrete method (NADM) and its improved version for suppressing numerical dispersion were developed recently. This new method is a refinement of two previous methods and further increases the ability of suppressing numerical dispersion for modeling acoustic wave propagation in 2D heterogeneous media, which uses acoustic wave displacement, particle velocity, and their gradients to restructure the acoustic wavefield via the truncated Taylor expansion and the high-order interpolation approximate method. For the method proposed, we investigate its implementation and compare it with the higher-order Lax-Wendroff correction (LWC) scheme, the original nearly analytic discrete method (NADM) and its im-proved version with regard to numerical dispersion, computational costs, and computer storage requirements. The numerical dispersion relations provided by the refined algorithm for 1D and 2D cases are analyzed, as well as the numerical results obtained by this method against the exact solution for the 2D acoustic case. Numerical results show that the refined method gives no visible numerical dispersion for very large spatial grid increments. It can simulate high-frequency wave propagation for a given grid interval and automatically suppress the numerical dispersion when the acoustic wave equation is discretized, when too few samples per wavelength are used, or when models have large velocity contrasts. Unlike the high-order LWC methods, our present method can save substantial computational costs and memory requirements because very large grid increments can be used. The refined method can be used for the simulation of large-scale wavefields.

Forward Semi-Lagrangian Advection with Mass Conservation and Positive Definiteness for Falling Hydrometeors

Abstract A semi-Lagrangian advection scheme is developed for falling hydrometeors in hopes of replacing the conventional Eulerian scheme that has been widely used in the cloud microphysics scheme of numerical atmospheric models. This semi-Lagrangian scheme uses a forward advection method to determine the advection path with or without iteration, and advected mass in a two-time-level algorithm with mass conservation. Monotonicity is considered in mass-conserving interpolation between Lagrangian grids and model Eulerian grids, thus making it a positive definite advection scheme. For mass-conserving interpolation between the two grid systems, the piecewise constant method (PCM), piecewise linear method (PLM), and piecewise parabolic method (PPM) are proposed. The falling velocity at the bottom cell edge is modified to avoid unphysical deformation by scanning from the top layer to the bottom of the model, which enables the use of a large time step with reasonable accuracy. The scheme is implemented and tested in the Weather Research and Forecasting (WRF) Single-Moment 3-Class Microphysics Scheme (WSM3). In a theoretical test bed with constant terminal velocity, the proposed semi-Lagrangian algorithm shows that the higher-order interpolation scheme produces less diffusive features at maximal precipitation. Results from another idealized test bed with mass-weighted terminal velocity demonstrate that the accuracy of the proposed scheme is still satisfactory even with a time step of 120 s when the mean terminal velocity averaged at the departure and arrival points is employed. A two-dimensional (2D) squall-line test using the WSM3 scheme shows that the control run with the Eulerian advection scheme and the semi-Lagrangian run with the PCM method reveal similar results, whereas behaviors using the PLM and PPM are similar with higher-resolution features, such as mammatus-like clouds.

High-order parabolic beam approximation for aero-optics

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

Improvement of conservative property and interface resolution of mesh-based two-phase flow simulation algorithms for splashing fluid behaviors

This study is concerned with improvement of numerical code called CRIMSON (Civa RefIned Multiphase SimulatiON), which has been developed to evaluate multi-phase flow behaviors based on the recent CFD (computational fluid dynamics) techniques. The CRIMSON employs a finite-volume method combined with the high order interpolation scheme, CIVA (cubic-interpolation with area/volume coordinates). The CRIMSON solves gas–liquid two phases by a unified scheme of CUP (combined unified procedure). The conventional CIVA method has two problems of losing interface sharpness in long-term calculation and non-conservativeness. In this study, the problems were solved by introduction of the ideas of the level set method and the phase field method. We verified our method by applying it to some popular benchmark problems of single bubble rising and collapse of water column problems.

A posteriori inclusion of parton density functions in NLO QCD final-state calculations at hadron colliders: the APPLGRID project

A method to facilitate the consistent inclusion of cross-section measurements based on complex final-states from HERA, TEVATRON and the LHC in proton parton density function (PDF) fits has been developed. This can be used to increase the sensitivity of LHC data to deviations from Standard Model predictions. The method stores perturbative coefficients of NLO QCD calculations of final-state observables measured in hadron colliders in look-up tables. This allows the a posteriori inclusion of parton density functions (PDFs), and of the strong coupling, as well as the a posteriori variation of the renormalisation and factorisation scales in cross-section calculations. The main novelties in comparison to original work on the subject are the use of higher-order interpolation of Lagrangian form, which substantially improves the trade-off between accuracy and memory use, and a CPU and computer memory optimised way to construct and store the look-up table using modern software tools. It is demonstrated that a sufficient accuracy on the cross-section calculation can be achieved with reasonably small look-up table size by using the examples of jet production and electro-weak boson (Z, W) production in proton-proton collisions at a center-of-mass energy of 14 TeV at the LHC. The use of this technique in PDF fitting is demonstrated in a PDF-fit to HERA data and simulated LHC jet cross-sections as well as in a study of the jet cross-section uncertainties at various centre-of-mass energies.

Constrained approximation of rational Bézier curves based on a matrix expression of its end points continuity condition

For high order interpolations at both end points of two rational Bézier curves, we introduce the concept of C ( v , u ) -continuity and give a matrix expression of a necessary and sufficient condition for satisfying it. Then we propose three new algorithms, in a unified approach, for the degree reduction of Bézier curves, approximating rational Bézier curves by Bézier curves and the degree reduction of rational Bézier curves respectively; all are in L 2 norm and C ( v , u ) -continuity is satisfied. The algorithms for the first and second problems can get the best approximation results, and for the third one, resorting to the steepest descent method in numerical optimization obtains a series of degree reduced curves iteratively with decreasing approximation errors. Compared to some well-known algorithms for the degree reduction of rational Bézier curves, such as the uniformizing weights algorithm, canceling the best linear common divisor algorithm and shifted Chebyshev polynomials algorithm, the new one presented here can give a better approximation error, do multiple degrees of reduction at a time and preserve high order interpolations at both end points.

A HIGH ACCURACY CONFORMAL METHOD FOR EVALUATING THE DISCONTINUOUS FOURIER TRANSFORM

A highly accurate, fast algorithm is proposed to evaluate the flnite Fourier transform of both continuous and discontinues functions. As the discretization is conformal to the function discontinuities, this method is called the conformal Fourier transform (CFT) method. It is applied to computational electromagnetics to calculate the Fourier transform of induced electric current densities in a volume integral equation. The spectral discrimination in the CFT method can be arbitrary and the spectral range can be as large as needed. As no discretization for the Fourier exponential kernel is needed, the CFT method is not restricted by the Nyquist sampling theorem, thus avoiding the aliasing distortions that exist in other traditional methods. The accuracy of the CFT method is greatly improved since the method is based on high order interpolation and the closed-form Fourier transforms for polynomials partly reduce the error due to discretization. Assuming Ns and N are the numbers

Development of Noncentered Wavenumber-Based Optimized Interpolation Schemes with Amplification Control for Overlapping Grids

A noncentered interpolation technique has been constructed to perform simulations using overlapping grids for complex geometries. High-order centered Lagrange polynomial interpolations and interpolations optimized in the Fourier space are first generalized to the noncentered case. These noncentered interpolations either generate significant dispersion errors or strongly amplify high-wavenumber components. Accordingly, a noncentered high-order wavenumber-based optimized interpolation method is developed with the addition of a nonlinear constraint for the control of the amplitude amplification induced by decentering. High-order piecewise polynomial regressions of the obtained interpolation coefficients are performed. The time stability of the method is investigated in the 1-D case when the interpolation method is used in conjunction with explicit high-order differencing, filtering schemes, as well as a 6-step Runge–Kutta time integration algorithm. A criterion is formulated to predict its stability as a function of the filtering strength and the Courant–Friedrichs–Lewy constant. Finally, 1-D convection simulations are presented to illustrate the stability and the accuracy of the developed noncentered interpolations.

A Rational Interpolation Scheme with Superpolynomial Rate of Convergence

The purpose of this study is to construct a high-order interpolation scheme for arbitrary scattered datasets. The resulting function approximation is an interpolation function when the dataset is exact, or a regression if measurement errors are present. We represent each datapoint with a Taylor series, and the approximation error as a combination of the derivatives of the target function. A weighted sum of the square of the coefficient of each derivative term in the approximation error is minimized to obtain the interpolation approximation. The resulting approximation function is a high-order rational function with no poles. When measurement errors are absent, the interpolation approximation converges to the target function faster than any polynomial rate of convergence.

Advection scheme with 3rd high-order spatial interpolation at the middle temporal level and its application to saltwater intrusion in the Changjiang Estuary

A non-oscillatory advection scheme with 3rd high-order spatial interpolation at the middle temporal level coupled with a TVD limiter (the 3rd HSIMT-TVD scheme) was developed with the following features: (1) It uses a parabola to characterize a substance’s local distribution at nΔ t; (2) it determines the cell interface value at ( n + 0.5)Δ t by one-step backward tracing to its position at nΔ t; and (3) it uses Sweby’s TVD limiter to constrain the derived cell interface value to eliminate numerical oscillation. The 3rd HSIMT-TVD scheme was derived based on Sweby’s original TVD formula directly, thus the symmetric assumption of a TVD limiter is not necessary. 1D/2D tests suggest that the present scheme performs quite robustly compared to some other well-established schemes. The wavenumber two freshwater eddy instability has been successfully simulated using the present scheme. The derived model was then used to simulate the saltwater intrusion into the Changjiang Estuary, and the in situ data collected in February 2003 and February–March 2007 were used for validation. The model results were quite close to the measured data, and the main features of saltwater intrusion, such as saltwater spilling over from the North Branch, the salt wedge, and the stratification introduced by the wedge-like saltwater intrusion and the plume-like extension of diluted water, were simulated correctly. Furthermore, some other schemes were also used for comparison, and their differences are discussed.

Multiscale computations for carbon nanotubes based on a hybrid QM/QC (quantum mechanical and quasicontinuum) approach

An effective multiscale computing scheme based on QM/QC (quantum mechanics/quasicontinuum) is applied for simulation of Carbon nanotubes (CNTS) mechanics. First, quasicontinuum simulation of deformations of curved crystalline structures is conducted to examine the fully nonlocal behavior of CNTs with the aid of high-order interpolation functions and the “cluster” concept, which facilitates accurate energy approximation for crystals. Next, a multiscale computing approach based on QM/QC hybridization is devised, and applied for simulation of CNT mechanics. The bonding configuration changes, e.g. bond breaking or creation, near defect sites are correctly represented with the QM/QC hybrid model. For studying electronic properties coupled with the mechanical deformation of CNTs, the change of the electrical properties from an initial semiconductor into metal under mechanical bending is investigated. Single-walled CNTs having various types of defects and subjected to uniaxial tension are considered for fracture. The theoretical strength of the CNTs in the presence of each defect is computed based on the QM/QC hybrid scheme, wherein the defect neighborhood is modeled as a QM zone for a first-principle-based calculation using density functional theory (DFT), and the remaining area as a QC zone. This multiscale computing approach greatly improves the accuracy in the prediction of the failure strains of CNTs over a purely molecular mechanical or quasicontinuum method.

A HIGH-ORDER FINITE ELEMENT FORMULATION FOR VIBRATION ANALYSIS OF BEAM-TYPE STRUCTURES

A high-order finite element model is presented to perform the vibration analysis of beams. The equations of motion are formulated by applying the principle of total potential energy in elastic dynamic system and the "set-in-right-position" rule for the construction of system matrices first proposed by the author. The primary advantage of the principle and rule lies in its simplicity and efficiency in solving the modeling problem of complex dynamic system. The requirement of strain continuity has certainly not being met at element interfaces with the use of conventional cubic Hermitian formulation. Hence, it is difficult to predict the dynamic responses of beams accurately. In order to overcome this problem, a beam element with simple higher-order interpolation function is chosen as the analysis model. Although the number of nodal degrees of freedom is increased herein, usually a coarse mesh will suffice. The present formulation is able to provide results of high accuracy with low computational effort. For the purpose of illustration, the dynamic characteristics analysis and dynamic response analysis are carried out on beam models. The solutions obtained for all the examples are in good agreement with the exact solutions found by fundamental theory of vibration.

A numerical source of small-scale number-density fluctuations in Eulerian–Lagrangian simulations of multiphase flows

Eulerian–Lagrangian simulations of multiphase flow are known to suffer from two errors that can introduce small-scale fluctuations in the number-density of the dispersed phase. These errors can be reduced by increasing the number of particles in the simulation. Here, we present results to demonstrate that a third error exists that can also generate small-scale number-density fluctuations. In contrast to the two known errors, this error cannot be lowered by increasing the number of particles. Analysis shows that this error is caused by spatial variation at the subgrid scale in the interpolation error of the fluid velocity to the particle location. If the particle velocity divergence is zero, the particle concentration error remains at the subgrid scale. However, if particles preferentially accumulate either due to their inertia or due to divergence of the underlying fluid-velocity field, this error manifests as number-density fluctuations on the grid scale. The only mechanism of reducing these errors is through higher-order accurate interpolation. By studying two model problems, estimates for the errors are derived. These estimates are shown to be quite accurate for simulations of shock and expansion waves interacting with particles.