High-order Numerical Methods Research Articles

Accurate, high-order representation of complex three-dimensional surfaces via Fourier continuation analysis

We present a new method for construction of high-order parametrizations of surfaces: starting from point clouds, the method we propose can be used to produce full surface parametrizations (by sets of local charts, each one representing a large surface patch – which, typically, contains thousands of the points in the original point-cloud) for complex surfaces of scientific and engineering relevance. The proposed approach accurately renders both smooth and non-smooth portions of a surface: it yields super-algebraically convergent Fourier series approximations to a given surface up to and including all points of geometric singularity, such as corners, edges, conical points, etc. In view of their C ∞ smoothness (except at true geometric singularities) and their properties of high-order approximation, the surfaces produced by this method are suitable for use in conjunction with high-order numerical methods for boundary value problems in domains with complex boundaries, including PDE solvers, integral equation solvers, etc. Our approach is based on a very simple concept: use of Fourier analysis to continue smooth portions of a piecewise smooth function into new functions which, defined on larger domains, are both smooth and periodic. The “continuation functions” arising from a function f converge super-algebraically to f in its domain of definition as discretizations are refined. We demonstrate the capabilities of the proposed approach for a number of surfaces of engineering relevance.

High-order numerical method for the nonlinear Helmholtz equation with material discontinuities in one space dimension

The nonlinear Helmholtz equation (NLH) models the propagation of electromagnetic waves in Kerr media, and describes a range of important phenomena in nonlinear optics and in other areas. In our previous work, we developed a fourth order method for its numerical solution that involved an iterative solver based on freezing the nonlinearity. The method enabled a direct simulation of nonlinear self-focusing in the nonparaxial regime, and a quantitative prediction of backscattering. However, our simulations showed that there is a threshold value for the magnitude of the nonlinearity, above which the iterations diverge. In this study, we numerically solve the one-dimensional NLH using a Newton-type nonlinear solver. Because the Kerr nonlinearity contains absolute values of the field, the NLH has to be recast as a system of two real equations in order to apply Newton’s method. Our numerical simulations show that Newton’s method converges rapidly and, in contradistinction with the iterations based on freezing the nonlinearity, enables computations for very high levels of nonlinearity. In addition, we introduce a novel compact finite-volume fourth order discretization for the NLH with material discontinuities. Our computations corroborate the design fourth order convergence of the method. The one-dimensional results of the current paper create a foundation for the analysis of multidimensional problems in the future.

High order numerical methods to a type of delta function integrals

We study second to fourth order numerical methods to a type of delta function integrals in one to three dimensions. These delta function integrals arise from recent efficient level set methods for computing the multivalued solutions of nonlinear PDEs. We show that the natural quadrature approach with usual discrete delta functions and support size formulas to the two dimensional delta function integrals suffer from nonconvergence. We then design high order numerical methods to this type of delta function integrals based on interpolation approach. Numerical examples are presented to verify the efficiency and accuracy of our methods.

High-order methods and numerical boundary conditions

In this paper we present high-order difference schemes for convection diffusion problems. When we apply high-order numerical methods to problems where physical boundary conditions are not periodic there is a need to choose adequate numerical boundary conditions in order to preserve the high-order accuracy. Next to the boundary we do not usually have enough discrete points to apply the high-order scheme and therefore at these nodes we must consider different approximations, named the numerical boundary conditions. The choice of numerical boundary conditions can influence the overall accuracy of the scheme and most of the times do influence the stability. Here, we discuss which orders of accuracy are reasonable to be considered at the numerical boundary conditions, such that we do not pay a high price in accuracy and stability.

NUMERICAL STUDY OF THE EFFECT OF INJECTION MODELING ON MIXING AND FLAME STABILITY IN SUPERSONIC COMBUSTION

The present study describes a numerical study on mixing and flame stability of hydrogen in supersonic airflow using injection modelling. 3 cases of injection model including 1 injector, 2 injectors and 3 injectors are considered. For these cases, numerical results are obtained by solving N-S equations, turbulence model and finite rate full chemistry of hydrogen-air using a high order numerical method. Temperature, flow field and mole fraction of all species are derived for all 3 cases. The results show a strong change in flame configuration for the 3 cases due to strong change in flow field.

A spectral finite volume transport scheme on the cubed-sphere

Advective processes are of central importance in many applications and their treatment is crucial in the numerical modelling of the transport of trace constituents in atmospheric models. High-order numerical methods offer the promise of accurately capturing these advective processes in atmospheric flows and have been shown to efficiently scale to large numbers of processors. In this paper, a conservative transport scheme based on the nodal high-order spectral finite volume method is developed for the cubed-sphere. A third-order explicit strong stability preserving scheme is employed for the time integration. The reconstruction procedure which we developed avoids the (expensive) calculation of the inverse of the reconstruction matrix. Flux-corrected transport algorithm is implemented to enforce monotonicity in the two-dimensional transport scheme. Two standard advection tests, a solid-body rotation and a deformational flow, were performed to evaluate the spectral finite volume method optionally combined with a flux-corrected transport scheme. Spectral accuracy in space is demonstrated with a linear wave equation.

Direct Numerical Simulation of a Non-Premixed Impinging Jet Flame

A non-premixed impinging jet flame at a Reynolds number 2000 and a nozzle-to-plate distance of two jet diameters was investigated using direct numerical simulation (DNS). Fully three-dimensional simulations were performed employing high-order numerical methods and high-fidelity boundary conditions to solve governing equations for variable-density flow and finite-rate Arrhenius chemistry. Both the instantaneous and time-averaged flow and heat transfer characteristics of the impinging flame were examined. Detailed analysis of the near-wall layer was conducted. Because of the relaminarization effect of the wall, the wall boundary layer of the impinging jet is very thin, that is, in the regime of viscous sublayer. It was found that the law-of-the-wall relations for nonisothermal flows in the literature need to be revisited. A reduced wall distance incorporating the fluid dynamic viscosity was proposed to be used in the law-of-the-wall relations for nonisothermal flows, which showed improved prediction over the law of the wall with the reduced wall distance defined in terms of fluid kinematic viscosity in the literature. Effects of external perturbation on the dynamic behavior of the impinging flame were found to be insignificant.

High-order numerical solution of the nonlinear Helmholtz equation with axial symmetry

The nonlinear Helmholtz (NLH) equation models the propagation of intense laser beams in a Kerr medium. The NLH takes into account the effects of nonparaxiality and backward scattering that are neglected in the more common nonlinear Schrödinger model. In [G. Fibich, S. Tsynkov, High-order two-way artificial boundary conditions for nonlinear wave propagation with backscattering, J. Comput. Phys., 171 (2001) 632–677] and [G. Fibich, S. Tsynkov, Numerical solution of the nonlinear Helmholtz equation using nonorthogonal expansions, J. Comput. Phys., 210 (2005) 183–224], a novel high-order numerical method for solving the NLH was introduced and implemented in the case of a two-dimensional Cartesian geometry. The NLH was solved iteratively, using the separation of variables and a special nonlocal two-way artificial boundary condition applied to the resulting decoupled linear systems. In the current paper, we propose a major improvement to the previous method. Instead of using LU decomposition after the separation of variables, we employ an efficient summation rule that evaluates convolution with the discrete Green's function. We also extend the method to a three-dimensional setting with cylindrical symmetry, under both Dirichlet and Sommerfeld-type transverse boundary conditions.

A staggered compact finite difference formulation for the compressible Navier–Stokes equations

In this paper, a compact high order (up to 12th order) numerical method to solve the compressible Navier–Stokes equations will be presented. A staggered arrangement of the variables has been used. It is shown that the method is not only very accurate but numerically also very stable even in the case that not all the energy containing scales in the flow are resolved. This in contrast to standard (collocated) compact finite difference methods. Some results for a turbulent non-reacting and a reacting jet with a Reynolds number of 10,000 and a Mach number of 0.5 are reported.

High-order solutions of three-dimensional rough-surface scattering problems at high frequencies. I: The scalar case

We present a new high-order numerical method for the solution of high-frequency scattering problems from rough surfaces in three dimensions. The method is based on the asymptotic solution of appropriate integral equations in the high-frequency regime, in a manner that bypasses the need to resolve the fields on the scale of the wavelength of radiation. Indeed, inspired by prior work in two dimensions, we seek a solution of the integral equation in the form of a slow modulation of the incoming radiation, and we choose a series expansion in inverse powers of the wavenumber to represent the unknown slowly varying envelope. As we show, this framework can be made to yield an efficiently computable recursion for the terms in the series to any arbitrary order. The resulting algorithms generally provide a very significant improvement over classical (e.g. Kirchhoff) approximations in both accuracy and applicability and they can, in fact, effectively produce results with full double-precision accuracy for configurations of practical interest and up to the resonance regime.

High-order shock-capturing methods for modeling dynamics of the solar atmosphere

We use one-dimensional high-order central shock-capturing numerical methods to study the response of various model solar atmospheres to forcing at the solar surface. The dynamics of the atmosphere is modeled with the Euler equations in a variable-sized flux tube in the presence of gravity. We study dynamics of the atmosphere suggestive of spicule formation and coronal oscillations. These studies are performed on observationally derived model atmospheres above the quiet sun and above sunspots. To perform these simulations, we provide a new extension of existing second- and third-order shock-capturing methods to irregular grids. We also solve the problem of numerically maintaining initial hydrostatic balance via the introduction of new variables in the model equations and a careful initialization mechanism. We find several striking results: all model atmospheres respond to a single impulsive perturbation with several strong shock waves consistent with the rebound-shock model. These shock waves lift material and the transition region well into the initial corona, and the sensitivity of this lift to the initial impulse depends nonlinearly on the details of the atmosphere model. We also reproduce an observed 3 min coronal oscillation above sunspots as well as 5 min oscillations above the quiet sun.

Axis switching in turbulent buoyant diffusion flames

Dynamics of buoyant diffusion flames from rectangular, square, and round fuel sources were investigated using direct numerical simulation (DNS). Fully three-dimensional simulations were performed employing high-order numerical methods and boundary conditions to solve governing equations for variable-density flow and finite-rate Arrhenius chemistry. Significant differences among the different cases were revealed in the vortex dynamics, entrainment rate, small-scale mixing, and consequently flame structures. Mixing and entrainment enhancement in non-circular flames in comparison with circular ones was explained using the Biot–Savart instability theory, which relates vortex dynamics to the local azimuthal curvature. An extension of the theory elucidated why rectangular flames entrain more efficiently and spread wider than square ones, although both configurations have corners. It also provided an explanation for the aspect ratio effects in the near field. In the far field, nonlinear effects were dominant and the general transport equations for vorticity were analyzed in detail. The corner effects and aspect ratio effects were shown to be augmented by the intricate interactions among vortex dynamics, combustion, and buoyancy through the various terms in the equations. The presence of corners in non-circular flames led to concentrated regions of fine-scale mixing and intense reactions centered around the corners. Moreover, the rectangular flames exhibited a different dynamic behavior from even the square one, by creating discrepancies in entrainment, mixing, and combustion between the minor and major axis directions. Increasing the aspect ratio exacerbated such directional discrepancies, and ultimately led to axis switching. It was the first time that axis switching was observed by DNS in a rectangular flame of aspect ratio 3, which raised further questions in combustion prediction and control. Finally, a unified explanation for corner and aspect ratio effects was given on the basis of the Biot–Savart instability theory and the vorticity transport equations.

Least squares spectral element method for 2D Maxwell equations in the frequency domain

We propose a high order numerical method for the first order Maxwell equations in the frequency domain, defined in media with arbitrary complex shape. Our approach is based on the combination of th...

Least squares spectral element method for 2D Maxwell equations in the frequency domain

We propose a high order numerical method for the first order Maxwell equations in the frequency domain, defined in media with arbitrary complex shape. Our approach is based on the combination of the least squares approach with the spectral element method. The former frees the solution from spurious modes, that can be found sometimes in classical finite element simulations. Many examples of such non-physical solutions exist in literature, and elimination of these spurious effects is a subject of great interest. Spectral elements are a numerical technique for solving partial differential equations which can be regarded as an extension of finite elements: they merge the flexibility of finite elements in dealing with complex geometries, and the better accuracy of spectral methods. Convergence to exact solution is improved by increasing (at run time) the polynomial degree, with no changes on the computational grid: this provides a significant advantage in respect to low order finite element, which necessarily have to resort to grid refinement. In the authors opinion this approach can be successfully used for the treatment of large scale electromagnetic problems or, alternatively, for applications where higher precision is required. We present a few numerical experiments which prove the capability of the method in object. Copyright © 2004 John Wiley & Sons, Ltd.

Numerical methods for stochastic differential equations.

Stochastic differential equations (SDE's) play an important role in physics but existing numerical methods for solving such equations are of low accuracy and poor stability. A general strategy for developing accurate and efficient schemes for solving stochastic equations is outlined here. High-order numerical methods are developed for the integration of stochastic differential equations with strong solutions. We demonstrate the accuracy of the resulting integration schemes by computing the errors in approximate solutions for SDE's which have known exact solutions.

High‐order numerical methods for one‐dimensional parabolic singularly perturbed problems with regular layers

AbstractIn this work we construct and analyze some finite difference schemes used to solve a class of time‐dependent one‐dimensional convection‐diffusion problems, which present only regular layers in their solution. We use the implicit Euler or the Crank‐Nicolson method to discretize the time variable and a HODIE finite difference scheme, defined on a piecewise uniform Shishkin mesh, to discretize the spatial variable. In both cases we prove that the numerical method is uniformly convergent with respect to the diffusion parameter, having order near two in space and order one or 3/2, depending on the method used, in time. We show some numerical examples which illustrate the theoretical results, in the case of using the Euler implicit method, and give better numerical behaviour than that predicted theoretically, showing order two in time and order N−2log2N in space, if the Crank‐Nicolson scheme is used to discretize the time variable. Finally, we construct a numerical algorithm by combining a third order A‐stable SDIRK with two stages and a third‐order HODIE difference scheme, showing its uniformly convergent behavior, reaching order three, up to a logarithmic factor. © 2004 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2005

Instabilities, entrainment and mixing in reacting plumes

Dynamics of buoyant reacting flows from rectangular, square and round sources were investigated using Direct Numerical Simulation (DNS). The general three-dimensional time-dependent governing equations for compressible flow and finite-rate Arrhenius chemistry were solved by high-order numerical methods. The study focused on the intricate couplings among instabilities, vortex dynamics, mixing, entrainment, turbulence and combustion through buoyancy. The main dynamic features of buoyant reacting flows, including the puffing phenomenon, were well reproduced. The instability behind puffing was identified as an intrinsic buoyancy instability that arose from density inhomogeneity in the presence of gravity. No external disturbances were required to trigger or sustain the puffing motions, indicating that the buoyancy instability was a global, absolute instability. The base configurations at the inlet had a significant influence on the subsequent vortex dynamics, entrainment, mixing and ultimately the combustion processes. Base configurations with corners, such as squares and rectangles, were associated with higher entrainment rates than those associated with circular jets/plumes, due to the Biot–Savart instability. A new mechanism based on an extension of the Biot–Savart instability clearly elucidated why rectangular jets/plumes entrain more than square ones, despite the fact that both types of configurations have corners. The aspect ratio effects arising from the mechanism could explain the higher level of vorticity and consequently higher entrainment rate along the major axis as compared with those along the minor axis. Axis switching, as may be deduced from the mechanism, was observed in the rectangular case of aspect ratio 3. The main terms in the vorticity transport equations were calculated for this case, and the coupling between vortex dynamics and combustion was scrutinized. The base configurations were seen to trigger different dynamic couplings among instabilities, vortex dynamics, large-scale entrainment, small-scale mixing and combustion. Such couplings were truly two-way, as exemplified by the modification of the energy frequency spectra by chemical heat release.

Large-eddy Simulation of Supersonic Boundary-layer Flow by a High-order Method

A high-order numerical method is described for performing large-eddy simulations (LESs) of supersonic flowfields. Spatial derivatives are represented by a compact stencil that is used in conjunction with a tenth-order non-dispersive filter. The scheme employs a time-implicit approximately-factored finite-difference algorithm, and applies Newton-like subiterations to achieve second-order temporal and sixth-order spatial accuracy. Details of the method are summarized and LESs are carried out for a spatially evolving supersonic turbulent boundary layer. Two different subgrid-scale models are incorporated in the simulations to account for the spatially under-resolved stresses and heat flux. Comparisons are made between the respective computations, as well as with available experimental data and with previous numerical results.

Large-Eddy Simulation of Supersonic Cavity Flowfields Including Flow Control

Large-eddy simulations of supersonic cavity flowfields are performed using a high-order numerical method. Spatial derivatives are represented by a fourth-order compact approximation that is used in conjunction with a sixth-order nondispersive filter. The scheme employs a time-implicit approximately factored finite difference algorithm, and applies Newton-like subiterations to achieve second-order temporal and fourth-order spatial accuracy. The Smagorinsky dynamic subgrid-scale model is incorporated in the simulations to account for the spatially underresolved stresses. Computations at a freestream Mach number of 1.19 are carried out for a rectangular cavity having a length-to-depth ratio of 5:1. The computational domain is described by 2.06×10 7 grid points and has been partitioned into 254 zones, which were distributed on individual processors of a massively parallel computing platform. Active flow control is applied through pulsed mass injection at a very high frequency, thereby suppressing resonant acoustic oscillatory modes

Higher-Order Numerical Methods for Transient Wave Equations

Solving efficiently the wave equations involved in modeling acoustic, elastic or electromagnetic wave propagation remains a challenge both for research and industry. To attack the problems coming from the propagative character of the solution, the author constructs higher-order numerical methods to reduce the size of the meshes, and consequently the time and space stepping, dramatically improving storage and computing times. This book surveys higher-order finite difference methods and develops various mass-lumped finite (also called spectral) element methods for the transient wave equations, and presents the most efficient methods, respecting both accuracy and stability for each sort of problem. A central role is played by the notion of the dispersion relation for analyzing the methods. The last chapter is devoted to unbounded domains which are modeled using perfectly matched layer (PML) techniques. Numerical examples are given.