Numerical Dispersion Properties Research Articles

Optimized Traditional and Leapfrog CDI-FDTD Schemes With Reduced Dispersion and Numerical Properties Analysis

In order to improve the numerical dispersion property of the finite-difference time-domain (FDTD) method further, two optimized low-dispersion schemes are herein proposed, which are based on the complying-divergence implicit (CDI) FDTD method and its one-step leapfrog method. These proposed approaches effectively reduce the numerical dispersion errors by intervening with anisotropic parameters while maintaining unconditional stability. The numerical stability and dispersion properties are analyzed in detail to illustrate their effectiveness. The factors, which impact the numerical dispersion including propagation angles, Courant–Friedrichs–Lewy (CFL) number, nonuniform multiscale, mesh resolution, and frequency are carefully examined. In addition, by using the dual-/triple-band multiple-input multiple-output (MIMO) antennas as numerical examples, the practicality of the optimization schemes is comprehensively validated by comparing the calculation time, memory resource, and relative error with those from the traditional and leapfrog CDI-FDTD methods. It is worth celebrating the superiority of the proposed schemes not only in terms of numerical dispersion, but also in calculation efficiency, solution accuracy, and floating-point operands.

A rectangular-grid lattice spring model for modeling elastic waves in Poisson’s solids

The lattice spring model (LSM) combined with the velocity Verlet algorithm is a newly developed scheme for modeling elastic wave propagation in solid media. Unlike conventional wave equations, LSM is established on the basis of micromechanics of the subsurface media, which enjoys better dynamic characteristics of elastic systems. We develop a rectangular-grid LSM scheme for elastic waves simulation in Poisson’s solids, and the direction-dependent elastic constants are deduced to keep the isotropy of the discrete grid. The stability condition and numerical dispersion properties of LSM are discussed and compared with other numerical methods. The 2D and 3D numerical simulations are carried out using the rectangular-grid LSM, as well as the second- and fourth-order accuracy staggered finite-difference method (FDM). Wavefields obtained by LSM are fairly similar with those by analytical solution and FDM, which demonstrates the correctness of the proposed scheme and its capability of modeling elastic wave propagation in heterogeneous media. Moreover, we perform plane P-wave simulation through a semi-infinite cavity model via LSM and FDM, the recorded wavefield snapshots indicate that our proposed rectangular-grid LSM obtains more reasonable wavefield details compared with those of FDM, especially in media with high compliance and structure complexity. Our main contribution lies in offering an alternative simulation scheme for modeling elastic wave propagation in media with some kinds of complexities, which conventional FDM may fail to simulate.

Higher Order LOD-FDTD Methods and Their Numerical Dispersion Properties

Higher order locally one-dimensional finite-difference time-domain (HO-LOD-FDTD) methods are proposed here to improve the accuracy of the four-step LOD-FDTD, five-step LOD-FDTD, and six-step LOD-FDTD methods. In the proposed HO-LOD-FDTD methods, higher order finite-difference approximations are used for all the three space derivatives; therefore, they give improved accuracy in all the three directions. In addition to HO-LOD-FDTD methods, hybrid-HO LOD-FDTD (HHO-LOD-FDTD) methods are also proposed here. In the HHO-LOD-FDTD methods, higher order finite-difference approximations are used for the desired direction(s) only. As a result, they give improved accuracy in the desired direction(s), and take less computational time than the HO-LOD-FDTD methods. Moreover, to validate the analytical results, simulation results of the fourth-order (FO) and hybrid-FO LOD-FDTD methods are also discussed here.

Unconditionally stable time-domain elastic wave simulations by the alternative direction implicit method

The stability of the conventional staggered-grid finite-difference time-domain (FDTD) method for elastic wave simulations is limited by the Courant condition and material heterogeneity. Its computational efficiency is significantly hampered when the mesh size is much smaller than a wavelength (for geometric modeling accuracy) and/or with a high impedance contrast. An unconditionally stable alternating direction implicit (ADI) method is proposed to overcome the conditional stability. It is based on additively splitting the Crank Nicholson (CN) operator into two sub-operators and subsequently solving the CN scheme in two sub-steps. In each sub-step, a tri-diagonal matrix is formed based on one of the sub-operators and the associated field variables are solved implicitly with O(N) computational complexity. The rest of the field variables are solved explicitly. Due to its semi-implicit nature, it can be proved that the ADI method is unconditionally stable regardless of the time step size. The numerical dispersion properties of the ADI method are also analyzed theoretically. Several numerical examples of acoustic wave scattering from elastic targets are further provided to demonstrate its accuracy and efficiency.

An analysis of a discontinuous spectral element method for elastic wave propagation in a heterogeneous material

The numerical dispersion and dissipation properties of a discontinuous spectral element method are investigated in the context of elastic waves in one dimensional periodic heterogeneous materials. Their frequency dependence and elastic band characteristics are studied. Dispersion relations representing both pass band and stop band structures are derived and used to assess the accuracy of the numerical results. A high-order discontinuous spectral Galerkin method is used to calculate the complex dispersion relations in heterogeneous materials. Floquet---Bloch theory is used to derive the elastic band structure. The accuracy of the dispersion relation is investigated with respect to the spectral polynomial orders for three different cases of materials. Numerical investigations illustrate a spectral convergence in numerical accuracy with respect to the polynomial order based on the elastic band structure and a discontinuous jump of the maximum resolvable frequency within the pass bands resulting in a step-like increase of it with respect to the polynomial order.

AN EFFICIENT ALGORITHM FOR THE NOVEL WEAKLY CONDITIONALLY STABLE FDTD METHOD

In this paper, we present an efficient formulation of the novel weakly conditionally stable finite-difference time-domain (NWCS-FDTD) method for the electromagnetic problems with very fine structures in one or two directions. The formulation is obtained by using only algebraic manipulation of the original method, and therefore the numerical stability and dispersion properties can be preserved. Moreover, due to its simpler right-hand sides of the updating equations, the proposed algorithm is more efficient than the existing WCS-FDTD methods, allowing a significant reduction in the cost of CPU time. Numerical experiments are finally given to verify the accuracy and efficiency of the proposed method.

Discussion of paper ‘Development of a family of unconditionally stable explicit direct integration algorithms with controllable numerical energy dissipation’ by Chinmoy Kolay and James M. Ricles, Earthquake Engineering and Structural Dynamics 2014; 43:1361–1380

SummaryIt seems that the explicit KR‐α method (KRM), which was developed by Kolay and Ricles, is promising for the step‐by‐step integration because it simultaneously integrates unconditional stability, explicit formulation, and numerical dissipation together. It was shown that KRM can inherit the numerical dispersion and energy dissipation properties of the generalized‐α method [1] for a linear elastic system, and it reduces to CR method (CRM), which was developed by Chen and Ricles [2] if ρ∞ = 1 is adopted, where ρ∞ is the spectral radius of the amplification matrix of KRM as the product of the natural frequency and the step size tends to infinity. However, two unusual properties were found for KRM and CRM, and they might limit their application to solve either linear elastic or nonlinear systems. One is the lack of capability to capture the structural nonlinearity, and the other is that it is unable to realistically reflect the dynamic loading. Copyright © 2014 John Wiley & Sons, Ltd.

Analysis and Design of a New Dispersion Relation Preserving Alternate Direction Bidiagonal Compact Scheme

In the present work bidiagonal schemes have been analyzed first by using global spectral analysis for wave propagation problem. The properties given by the numerical amplification factor, dispersion and phase error are presented for these schemes to study effects on error propagation for 1D convection equation. Analysis of bidiagonal schemes helps in designing a new alternate direction bidiagonal (ADB) spatial scheme with improved numerical dispersion and dissipation properties. Design of the ADB scheme is facilitated by the prefactorization technique given in Hixon and Turkel (J Comput Phys 158:51---70, 2000) and is obtained here by minimizing the error in computing 1D convection equation using the correct stability and error analysis given in Sengupta et al. (J Comput Phys 226:1211---1218, 2007). The ADB scheme for spatial discretization is coupled with four-stage, fourth order accurate Runge---Kutta time marching scheme and optimized for multiple objective functions by classical tools to obtain a new space-time discretization scheme with desired numerical amplification and dispersion properties up to moderate wavenumbers, preserving physical dispersion relation almost exactly. The resultant optimal DRP scheme is fourth order accurate, both in space and time. We have compared this scheme with an accurate compact scheme in reproducing exact solution of 1D convection equation. Two benchmark problems from computational acoustics are shown for further verification of the accuracy of solution obtained by the present method. To demonstrate the practical utility of the method, we have presented the solution of zero pressure gradient flow past semi-infinite flat plate from the receptivity stage to fully developed inhomogeneous turbulent state without any model and compared with published results in Sengupta et al. (Phys Rev E 85:026308, 2012).

A Multiresolution Time-Domain Method for LEMP Calculation and Comparison With FDTD

A multiresolution time-domain (MRTD) method is proposed for cloud-to-ground lightning electromagnetic pulse (LEMP) calculation. Scaling function of Daubechies wavelet is chosen for fields expansion, and iterative formulas of MRTD method are derivated in 2-D cylindrical coordinates. First, three double-exponential lightning channel-base current with different rise times have been used for evaluating performance of the finite-difference time-domain (FDTD) method. And then, LEMP of different distances of the lightning channel were calculated with the MRTD method. We examined the MRTD method by comparing the results obtained with the FDTD method. Numerical simulation results show that the MRTD method and FDTD method are consistent with each other for first return stroke, even with a larger grid division. However, for subsequent return stroke, oscillations will be observed in the late response of the waveforms calculated by the FDTD method. Comparing with the FDTD method, the MRTD method suffers less oscillations with the same grid division, and can significantly save both computational memory and CPU time with large gird mesh, as well as good numerical dispersion property.

Hexagonal vs. rectilinear grids for explicit finite difference schemes for the two-dimensional wave equation

Finite difference schemes for the 2-D wave equation operating on hexagonal grids and the accompanying numerical dispersion properties have received little attention in comparison to schemes operating on rectilinear grids. This paper considers the hexagonal tiling of the wavenumber plane in order to show that the hexagonal grid is a more natural choice to emulate the isotropy of the Laplacian operator and the wave equation. Performance of the 7-pt scheme on a hexagonal grid is better than previously reported as long as the correct stability limit and tiling of the wavenumber plane are taken into account. Numerical dispersion is analyzed as a function of temporal frequency to demonstrate directional cutoff frequencies. A comparison to 9-pt compact explicit schemes on rectilinear grids is presented using metrics relevant to acoustical simulation. It is shown that the 7-pt hexagonal scheme has better computational efficiency than parameterized 9-pt compact explicit rectilinear schemes and that the error remains isotropic to fourth-order. Simulation results are presented.

Analysis of a mixed finite‐element pair proposed for an atmospheric dynamical core

AbstractWe present a numerical dispersion analysis for the P2 − P1DG finite‐element pair applied to the linear shallow‐water equations in one dimension. The aim is to provide insight into the numerical dispersion properties of the RT1 and BDFM1 finite‐element pairs in two dimensions, which have recently been proposed for horizontal discretisations of atmospheric dynamical cores with quasi‐uniform grids. This is achieved via analysis of a one‐dimensional RT1 element. Whilst these finite‐element pairs have been shown to have many desirable properties that extend properties of the C grid to non‐orthogonal quadrilateral and triangular grids, including stationary geostrophic modes on the f plane, and a 2:1 ratio of velocity to pressure degrees of freedom (a necessary condition for the absence of spurious mode branches), it is also important to have appropriately physical numerical wave propagation. In the absence of Coriolis force, we compute the group velocity for P2 − P1DG. We find that, as well as dropping to zero at the grid‐scale, which also occurs for the C‐grid finite‐difference method, the group velocity also drops to zero in a narrow band around kh = π which corresponds to eigenmodes with a wavelength close to two element widths. This is a potential problem because it increases the amount of wavenumber space that needs to be filtered. In this one‐dimensional case, we find that this particular issue can be removed by a small modification of the equations, namely partially lumping the mass matrix, in such a way that the other favourable properties of the scheme are not affected. We discuss various symmetric and asymmetric modifications of the mass matrix, and show that both such modifications preserve energy conservation (having modified the definition of discrete kinetic energy). Finally we illustrate our findings with numerical experiments, and discuss the potential to extend this modification to two‐dimensional schemes.

Dispersion analysis of the spectral element method

AbstractThe spectral element method (SEM) has (with exact time integration) the desirable attribute of locally and globally conserving mass, energy and potential vorticity. It also scales well on massively parallel computers. Another desirable attribute of a numerical method for an atmospheric dynamical core is that it should have good numerical dispersion properties in order to accurately represent wave propagation and adjustment processes. Application of the SEM to the one‐way wave equation is analysed to provide insight into its dispersion properties as a function of spectral order. For the lowest‐order spectral truncation (linear) the SEM discretisation is formally equivalent to centred second‐order finite differences on an Arakawa A grid. It consequently shares its poor dispersion properties, including energy propagation in the wrong direction for the short‐wavelength half of the spectrum. Increasing the spectral truncation of the SEM to quadratic improves its dispersion properties for the long‐wavelength part of the spectrum, but the problem of energy propagation in the wrong direction for the short‐wavelength part remains. Further increasing the order of the spectral truncation not only fails to address the poor energy propagation at small scales, but also introduces new problems, including gaps in the spectrum of frequencies that can be represented, and localisation of eigenmode structures near element boundaries. Numerical integrations confirm that these SEM dispersion properties lead to reversed group velocities and to grid imprinting at spectral element boundaries. Copyright © 2012 Royal Meteorological Society and British Crown Copyright, the Met Office

Coupled Discontinuous Galerkin/Finite Difference Solver on Hybrid Meshes for Computational Aeroacoustics

T HE field of direct simulation of acoustic waves propagation in the presence of a fluid in motion holds high standards of requirements. This is especially the case when it comes to industrial computational aeroacoustics (CAA), where methods must be able to copewith realistic, complex geometries while involving a high order of accuracy and yet result in affordable solvers in terms of computational resources. In that context, different computational methods have emerged over the years. Resulting from various theoretical backgrounds and approaches, they naturally hold specific characteristics, advantages, and drawbacks. Among them, discontinuous Galerkin (DG)methods [1–4] and finite difference (FD)methods [5– 7] are widely spread. The main characteristics of DG and FD methods are roughly recalled in Table 1. DGmethods arewell adapted to take into account complex geometries as they can deal with unstructured meshes [4,8]. Furthermore, their formulation naturally allows local-order refinement and is well suited for parallel computing. On the other hand, FD methods preferably run on structured meshes and are particularly efficient on Cartesian grids. They are easier to implement and less demanding in terms of CPU resources yet show good numerical dispersion and dissipation properties. Both these methods have been extensively studied, successfully implemented, and applied to that range of problems in industrial configurations. Based on this, a new family of questions have been raised. In particular, the possibility of coupling solvers of different natures in order to locally take advantage of the qualities of each method has already been studied byUtzmann et al. [9–11] in the field of CAA. We will also focus on hybridization techniques of DG/FD schemes based on a computational domain-decomposition approach. A typical motivation behind this is to be able to approximate the solution in the close neighborhood of complex obstacles on an unstructured DG mesh and compute the rest of the field on a Cartesian FD grid in order to alleviate computational time and resources. The purpose of this study is to investigate natural questions raised by such a coupling in order to gather informations on the behavior of the resulting hybrid solver in precision and stability. In this paper, we introduce two two-dimensional (2-D) DG timedomain (DGTD)/FD time-domain (FDTD) hybridization algorithms in the context of the resolution of the linearized Euler equations (LEEs) in two different kinds of geometrical configurations. Our results mainly focus on validation aspects, driven by numerical experiments that were conducted on academic test cases, as theoretical results on the interaction of these two schemes presently seem out of reach in a general case. We also present the behavior of our hybrid solver in the context of an acoustic benchmark problem that consists in the diffraction of an acoustic source by complex-geometries obstacles. The paper is organized as follows. In Sec. II, we present our physical modeling and both the FDTD andDGTD schemes that are used in this study.We also recall their compared properties, advantages, and drawbacks. In Sec. III, after presenting the numerical and theoretical issues raised by the design of a DG/FD coupled solver, we focus on our hybridization strategy and its numerical validation. In Sec. IV, we present an application of our hybrid solver to a complex-geometry acoustic test case, which was proposed at the Fourth CAA Workshop on Benchmark Problems [12].

Stability and dispersion analysis of the finite-difference time-domain algorithms in modelling the earth–ionosphere system

An analysis of the numerical stability and dispersion properties of the latitude–longitude and geodesic finite-difference time-domain (FDTD) algorithms in modelling the earth–ionosphere system is presented. To the best of the authors' knowledge, it is the first time the properties are analysed via rigorous theoretical derivations. Furthermore, both the analytical and experimental results are compared for the two algorithms. It is found that the geodesic FDTD algorithm has a looser stability limitation than that of the latitude–longitude FDTD algorithm. After the theoretical derivation and the numerical analysis, the global electromagnetic phenomenon was clearly exhibited using the two different algorithms, revealing the correctness of the theoretical analysis and the experiments. Additionally, the potential applications of the two algorithms are briefly proposed.

Inter-phase interaction in a turbulent, vertical channel flow laden with heavy particles. Part I: Numerical methods and particle dispersion properties

We investigated the interaction between two phases by numerical simulation of a turbulent channel flow laden with heavy particles. We calculated the fluid phase turbulence by a large eddy simulation with a point-force two-way coupling model and tracked particles by the Lagrangian trajectory method. Six cases with different particle volume concentrations for three kinds of particles were simulated. The particle dispersion properties on the inter-phase interaction were presented in this first part. The simulations provided comprehensive instantaneous two-phase turbulent flow fields. The particle streaky and cloudy preferential distributions, which weaken with increased particle inertia, are suppressed by particle modulation of turbulence. Modulation of turbulence is stronger near the wall than in other regions, because particle volume concentration is higher close to the wall and modulation increases with local particle concentration. The streaky particle distribution structures near the wall disappear earlier than the outer layer cloudy structures as the particle mass loading increases. Particle-phase instantaneous velocity scatter distributions, characterized by a Y-shape changing to a Λ-shape from the viscous sub-layer to the central plane, are altered by particle modulation of turbulence.

A Nonorthogonal ADI-FDTD Algorithm for Solving Two Dimensional Scattering Problems

In this paper, an alternating-direction implicit (ADI) scheme is applied to the finite-difference time-domain (FDTD) method for solving electromagnetic scattering problems in a generalized coordinate system. A formulation for two dimensional problems is presented and its numerical dispersion and stability property are discussed. In our generalized approach, the nonorthogonal grid is used to model the complex region of a scatterer only, whereas the standard FDTD lattice is used for the remaining regions. As a result, accurate griddings with a simple algorithm can be obtained using the new scheme, and the complexity of the algorithm is minimal. The perfectly matched layer (PML) is used to truncate the boundary. To illustrate the theory, a sinusoidal plane wave and a Gaussian pulse that propagates through a space modeled by locally nonorthogonal grids are used, with the stability of the code examined. The radar cross section of a perfectly conducting cylinder with a thin coating, a large curvature, and/or a sharp edge is calculated using the proposed method, and the result is compared with those using other conventional FDTD methods. It is found that the proposed algorithm is much more efficient than its FDTD counterpart when a complex object is analyzed.

Symplectic and multisymplectic Lobatto methods for the “good” Boussinesq equation

In this paper, we construct second order symplectic and multisymplectic integrators for the “good” Boussineq equation using the two-stage Lobatto IIIA-IIIB partitioned Runge–Kutta method, which yield an explicit scheme and is equivalent to the classical central difference approximation to the second order spatial derivative. Numerical dispersion properties and the stability of both integrators are investigated. Numerical results for different solitary wave solutions confirm the excellent long time behavior of symplectic and multisymplectic integrators by preserving local and global energy and momentum.

Determination of wave propagation properties for multilayered dielectric waveguide by semivectorial beam propagation method with Douglas scheme

AbstractAn improved parabolic‐equation‐based numerical method with Douglas scheme is developed for modeling the new 3D multilayered dielectric transmission waveguides. Unlike the traditional one, such waveguides are constructed by putting the multilayered dielectric cladding with step indices property on a silicon substrate. The propagation properties of widely used polarization‐maintaining waveguides and antiresonant reflecting dielectric waveguide are investigated by this method as examples. The discussion on the numerical dispersion and radiation loss property is also included to validate this improved beam propagation method (or BPM). © 2008 Wiley Periodicals, Inc. Microwave Opt Technol Lett 50: 2521–2524, 2008; Published online in Wiley InterScience (www.interscience.wiley.com). DOI 10.1002/mop.23725

Analysis of the numerical dispersion of waves in saturated poroelastic media

This work presents the analysis of the numerical dispersion properties of a finite element method used to approximate the solution of the equations of motion for a fluid saturated porous elastic solid (a Biot medium) in the two-dimensional case and the low frequency regime. The finite element method employed comprises a nonconforming rectangular element for the approximation of each component of the displacement vector in the solid phase, and the Raviart–Thomas–Nedelec mixed finite element space of zero order for the fluid phase. The study is carried out by constructing and analyzing analytic and numerical dimensionless dispersion relations, and by evaluating derived quantities such as dimensionless phase and group velocities and dimensionless attenuation for the three type of waves predicted by Biot’s equations of motion as a measure of the numerical distortion. It is observed that the finite element procedure introduces both numerical dispersion and anisotropy in all three waves in a similar fashion; however, only the slow wave displayed a significant frequency dependent behavior. It was also observed that the loss of accuracy is more important for the dimensionless attenuation than for the dimensionless group or phase velocities. The analysis presented yields lower bounds for the number of points per wavelength needed to reach a desired accuracy in the dimensionless phase and group velocities and attenuation coefficients in Biot media.

Direct Numerical Aeroacoustics for low Mach numbers using one Unstructured Grid

Using an unstructured grid, for which the mesh size is adapted to the local flow and acoustic solution, one can reduce the number of unknowns needed for a direct aeroacoustic simulation. The compressible Navier-Stokes equations for Mach numbers 0 <M <1 are considered. However, the unstructured grid together with the increasing stiffness of the problem for M → 0 reveals new numerical challenges. One is, to find an optimal time step size, that will reveal an equal numerical error distribution on the computational domain. Therefore, CFL-numbers greater than 1 are required, so that an implicit approach is chosen to overcome this complexity. The multigrid method used to solve the resulting algebraic system of equations is robust against the Mach number. Hence, the convergence rate does not deteriorate for M → 0. The accuracy of the time and space discretisation is investigated with respect to their numerical dispersion and diffusion properties. Furthermore, the order of the space discretisation is increased with a τ -extrapolation method, which is embedded into the multigrid procedure. The accuracy of the method and the reduction of unknowns is demonstrated by numerical simulations.