High-order Numerical Methods Research Articles

High order numerical methods for Vlasov-Poisson models of plasma sheaths

This article is a report of the CEMRACS 2022 project, called HIVLASHEA, standing for ”High order methods for Vlasov-Poisson models for sheaths”. A two-species Vlasov-Poisson model is described together with some numerical simulations, permitting to exhibit the formation of a plasma sheath. The numerical simulations are performed with two different methods: a first order classical finite difference (FD) scheme and a high order semi-Lagrangian (SL) scheme with Strang splitting; for the latter one, the implementation of (non-periodic) boundary conditions is discussed and different solvers for the electric field are compared. The codes are first evaluated on a one-species case, where an analytical solution is known. For the two-species case, cross-comparisons are performed and the numerical convergence of the SL method is investigated. The use of high order method is emphasized in this context. The (parallel) SL code turns out to be a robust tool to provide reliable numerical approximations of sheath stationary solution with realistic physical parameters.

A new $ \alpha $-robust nonlinear numerical algorithm for the time fractional nonlinear KdV equation

<abstract><p>In this work, we consider an $ \alpha $-robust high-order numerical method for the time fractional nonlinear Korteweg-de Vries (KdV) equation. The time fractional derivatives are discretized by the L1 formula based on the graded meshes. For the spatial derivative, the nonlinear operator is defined to approximate the $ uu_x $, and two coupling equations are obtained by processing the $ u_{xxx} $ with the order reduction method. Finally, the nonlinear difference schemes with order ($ 2-\alpha $) in time and order $ 2 $ precision in space are obtained. This means that we can get a higher precision solution and improve the computational efficiency. The existence and uniqueness of numerical solutions for the proposed nonlinear difference scheme are proved theoretically. It is worth noting the unconditional stability and $ \alpha $-robust stability are also derived. Moreover, the optimal convergence result in the $ L_2 $ norms is attained. Finally, two numerical examples are given, which is consistent with the theoretical analysis.</p></abstract>

A reduced-order Jacobi spectral collocation method for solving the space-fractional FitzHugh–Nagumo models with application in myocardium

PurposeThe space fractional PDEs (SFPDEs) play an important role in the fractional calculus field. Proposing a high-order, stable and flexible numerical procedure for solving SFPDEs is the main aim of most researchers. This paper devotes to developing a novel spectral algorithm to solve the FitzHugh–Nagumo models with space fractional derivatives.Design/methodology/approachThe fractional derivative is defined based upon the Riesz derivative. First, a second-order finite difference formulation is used to approximate the time derivative. Then, the Jacobi spectral collocation method is employed to discrete the spatial variables. On the other hand, authors assume that the approximate solution is a linear combination of special polynomials which are obtained from the Jacobi polynomials, and also there exists Riesz fractional derivative based on the Jacobi polynomials. Also, a reduced order plan, such as proper orthogonal decomposition (POD) method, has been utilized.FindingsA fast high-order numerical method to decrease the elapsed CPU time has been constructed for solving systems of space fractional PDEs.Originality/valueThe spectral collocation method is combined with the POD idea to solve the system of space-fractional PDEs. The numerical results are acceptable and efficient for the main mathematical model.

A robust higher-order fitted mesh numerical method for solving singularly perturbed parabolic reaction-diffusion problems

In this study, a robust higher-order numerical method for solving singularly perturbed parabolic reaction-diffusion problems is presented. The Crank-Nicolson method is applied to discretize the time derivative on a uniform mesh. On a Shishkin mesh, the space derivative is discretized using a hybrid numerical method that combines the cubic spline in tension method for the boundary layer regions with the central difference method for the outer region. Theoretically, we proved that the proposed hybrid numerical method is second-order in the outer region and fourth-order in the boundary layer regions in the space direction. As a result of this, the proposed method produces an almost second-order rate of convergence in the time domain and a higher-order rate of convergence in the space domain. The newly developed method is numerically demonstrated to be uniformly convergent at higher-order, independent of the perturbation parameter. Three numerical examples are computed to support the theoretical results.

A two-phase fluid model for epidemic flow.

We propose a new mathematical and computational modeling framework that incorporates fluid dynamics to study the spatial spread of infectious diseases. We model the susceptible and infected populations as two inviscid fluids which interact with each other. Their motion at the macroscopic level characterizes the progression and spread of the epidemic. To implement the two-phase flow model, we employ high-order numerical methods from computational fluid dynamics. We apply this model to simulate the COVID-19 outbreaks in the city of Wuhan in China and the state of Tennessee in the US. Our modeling and simulation framework allows us to conduct a detailed investigation into the complex spatiotemporal dynamics related to the transmission and spread of COVID-19.

Supersonic Base Flow With and Without Propulsive Jets: A Brief Review

Supersonic base flow is encountered in many aerospace applications related to satellite launch vehicle, missiles and projectiles. The origin and evolution of large scale structures which govern the dynamics of the base flow is not properly understood. The presence of central propulsive jet further complicates the base flow characteristics. Different semi-empirical, experimental and numerical studies carried out in supersonic base flow is briefly reviewed. Various issues and challenges faced in CFD simulation of base flow problems are described. Apart from resolution of numerics, proper choice of turbulence model remains the key aspect in the Reynolds Averaged Navier Stokes (RANS) simulations. Prediction of axial distribution of centerline velocities and recirculation bubble size remain problematic in RANS simulation. Detached Eddy Simulation (DES) results show dependence on model constants. Large Eddy Simulation (LES) studies with proper filter and Sub Grid Scale (SGS) model is promising, but its applications are limited to canonical problems due to very high computational requirements. Direct Numerical Simulation (DNS) studies reported in the literature consider much lower Reynolds number and hence could not be compared with the experimental data. It is suggested to use RANS and DES methods with anchored turbulence model from well-designed experiments for design data generation. Efforts should be made to use higher order numerical methods (LES and DNS) to explore the origin and evolution of large scale structures observed in compressible base flow.

The dissipative Generalized Hydrodynamic equations and their numerical solution

“Generalized Hydrodynamics” (GHD) stands for a model that describes one-dimensional integrable systems in quantum physics, such as ultra-cold atoms or spin chains. Mathematically, GHD corresponds to nonlinear equations of kinetic type, where the main unknown, a statistical distribution function f(t,z,θ), lives in a phase space which is constituted by a one-dimensional position variable z, and a one-dimensional “kinetic” variable θ, actually a wave-vector, called “rapidity”. Two key features of GHD equations are first a non-local and nonlinear coupling in the advection term, and second an infinite set of conserved quantities, which prevent the system from thermalizing. To go beyond this, we consider the dissipative GHD equations, which are obtained by supplementing the right-hand side of the GHD equations with a non-local and nonlinear diffusion operator or a Boltzmann-type collision integral. In this paper, we deal with new high-order numerical methods to efficiently solve these kinetic equations. In particular, we devise novel backward semi-Lagrangian methods for solving the advective part (the so-called Vlasov equation) by using a high-order time-Taylor series expansion for the advection fields, whose successive time derivatives are obtained by a recursive procedure. This high-order temporal approximation of the advection fields is used to design new implicit/explicit Runge–Kutta semi-Lagrangian methods, which are compared to Adams–Moulton semi-Lagrangian schemes. For solving the source terms, constituted by the diffusion and collision operators, we use and compare different numerical methods of the literature.

High-Order Spectral Method of Density Estimation for Stochastic Differential Equation Driven by Multivariate Gaussian Random Variables

There are some previous works on designing efficient and high-order numerical methods of density estimation for stochastic partial differential equation (SPDE) driven by multivariate Gaussian random variables. They mostly focus on proposing numerical methods of density estimation for SPDE with independent random variables and rarely research density estimation for SPDE is driven by multivariate Gaussian random variables. In this paper, we propose a high-order algorithm of gPC-based density estimation where SPDE driven by multivariate Gaussian random variables. Our main techniques are (1) we build a new multivariate orthogonal basis by adopting the Gauss–Schmidt orthogonalization; (2) with the newly constructed orthogonal basis in hand, we first assume the unknown function in the SPDE has the stochastic general polynomial chaos (gPC) expansion, second implement the stochastic gPC expansion for the SPDE in the multivariate Gaussian measure space, and third we obtain and numerical calculation deterministic differential equations for the coefficients of the expansion; (3) we used high-order algorithm of gPC-based for density estimation and moment estimation. We apply the newly proposed numerical method to a known random function, stochastic 1D wave equation, and stochastic 2D Schnakenberg model, respectively. All the presented stochastic equations are driven by bivariate Gaussian random variables. The efficiency is compared with the Monte-Carlo method based on the known random function.

High-order supplementary variable methods for thermodynamically consistent partial differential equations

In recent years, there has been growing interest in developing numerical methods to retain consistency between the partial differential equation system and its deduced equations (e.g., the thermodynamically consistent model and its energy law) after numerical approximations. A supplementary variable method (SVM) of second-order was recently introduced in [1,2] to preserve thermodynamic consistency for thermodynamically consistent models. In this paper, we present a new class of high-order supplementary variable methods. These methods consist of two parts: firstly, the original model is discretized using a high-order prediction method, and secondly, the SVM reformulated system is approximated by a high-order correction method. The newly developed schemes are rigorously proven to preserve high-order accuracy with increased stability through the use of a stabilizer. Additionally, the schemes inherently maintain thermodynamic consistency for thermodynamically consistent models, retaining the original energy law at the discrete level. To assess the high-order supplementary variable methods, we apply them to two thermodynamically consistent partial differential equation systems: the functionalized Cahn–Hilliard model and the ternary Cahn–Hilliard model. We then provide several benchmark examples to demonstrate the accuracy, stability, and efficiency of the newly developed high-order numerical methods.

Multi-dimensional summation-by-parts operators for general function spaces: Theory and construction

Summation-by-parts (SBP) operators allow us to systematically develop energy-stable and high-order accurate numerical methods for time-dependent differential equations. Until recently, the main idea behind existing SBP operators was that polynomials can accurately approximate the solution, and SBP operators should thus be exact for them. However, polynomials do not provide the best approximation for some problems, with other approximation spaces being more appropriate. We recently addressed this issue and developed a theory for one-dimensional SBP operators based on general function spaces, coined function-space SBP (FSBP) operators. In this paper, we extend the theory of FSBP operators to multiple dimensions. We focus on their existence, connection to quadratures, construction, and mimetic properties. A more exhaustive numerical demonstration of multi-dimensional FSBP (MFSBP) operators and their application will be provided in future works. Similar to the one-dimensional case, we demonstrate that most of the established results for polynomial-based multi-dimensional SBP (MSBP) operators carry over to the more general class of MFSBP operators. Our findings imply that the concept of SBP operators can be applied to a significantly larger class of methods than is currently done. This can increase the accuracy of the numerical solutions and/or provide stability to the methods.

Convergence of Collocation Methods for One Class of Impulsive Delay Differential Equations

This paper is concerned with collocation methods for one class of impulsive delay differential equations (IDDEs). Some results for the convergence, global superconvergence and local superconvergence of collocation methods are given. We choose a suitable piecewise continuous collocation space to obtain high-order numerical methods. Some illustrative examples are given to verify the theoretical results.

A high order numerical method for the variable order time-fractional reaction-subdiffusion equation

In this paper, we present a high order new numerical approximation for variable order Caputo fractional derivative of order 0<α(x,t)<1, by using the idea of interpolation. Then, by using this approximation, a numerical scheme is presented by using finite difference approach for variable order time fractional reaction-subdiffusion equation (VO-TFRSDE). The unconditionally stability of the numerical scheme is examined theoretically. The scheme is implemented on the two test problems. The numerical results are highly accurate with higher order of convergence.

Adaptive Model Reduction of High-Order Solutions of Compressible Flows via Optimal Transport

The solution of conservation laws with parametrised shock waves presents challenges for both high-order numerical methods and model reduction techniques. We introduce an r-adaptivity scheme based on optimal transport and apply it to develop reduced order models for compressible flows. The optimal transport theory allows us to compute high-order r-adaptive meshes from a starting reference mesh by solving the Monge–Ampère equation. A high-order discretization of the conservation laws enables high-order solutions to be computed on the resulting r-adaptive meshes. Furthermore, the Monge–Ampère solutions contain mappings that are used to reduce the spatial locality of the resulting solutions and make them more amenable to model reduction. We use a non-intrusive model reduction method to construct reduced order models of both the mesh and the solution. The procedure is demonstrated on three supersonic and hypersonic test cases, with the hybridisable discontinuous Galerkin method being used as the full order model.

Global instability phenomenon as a physical mechanism controlling dynamics of a nitrogen-diluted hydrogen flame

We investigate the dynamics of counter-current reacting flow under varying strengths of the counterflow and report the occurrence of global instability phenomena. This clarifies the role of global instability in the transition between attached and lifted flames. The research is performed with the help of the Large Eddy Simulation method using a high-order numerical method. The combustion process is modelled by applying detailed chemical kinetics for hydrogen combustion with chemical reaction terms computed from the resolved scale. The flow configuration consists of a central jet of nitrogen-diluted hydrogen fuel surrounded by an annular nozzle with a suction slot enforcing counterflow in the vicinity of the fuel stream. Suction strength is controlled by the ratio between the velocity in the suction slot (Usuc) and the velocity of the fuel jet (Uj) (I=Usuc/Uj). The Reynolds number based on the fuel parameters and the central nozzle diameter is Re=1600. Two values of the suction strength are considered (I=0.1,0.2) as well as the case without suction (I=0.0). A hot oxidizer (air) provided in the region outside the co-axial nozzle causes fuel auto-ignition far from the injection system after which the flame propagates downstream. It is shown that increasing the suction strength postpones the ignition and shifts its axial location downstream. Depending on the I parameter the flame stabilizes as attached (I=0.0,0.1) or lifted (I=0.2). It is shown that the flames in the cases with I=0.0,0.1 are very similar to each other appearing as laminar non-premixed flames both in the nozzle vicinity and far downstream. The flame in case I=0.2 is turbulent and much more dynamic. Its lifting is the result of global instability appearing at increased suction in the region of counter-current flow. It is shown that the global instability induces strong toroidal vortices that create a premixed reaction zone at the flame base. This flow structuring prevents upstream flame propagation and intensifies the mixing process ensuring almost complete fuel burning in a short distance from the nozzle.

An energy-conservative DG-FEM approach for solid–liquid phase change

We present a discontinuous Galerkin method for melting/solidification problems based on the “linearized enthalpy approach,” which is derived from the conservative form of the energy transport equation and does not depend on the use of a so-called mushy zone. We use the symmetric interior penalty method and the Lax–Friedrichs flux to discretize diffusive and convective terms, respectively. Time is discretized with a second-order implicit backward differentiation formula, and two outer iterations with second-order extrapolation predictors are used for the coupling of the momentum and energy. The numerical method was validated with three different benchmark cases, i.e., the one-dimensional Stefan problem, octadecane melting in a square cavity and gallium melting in a rectangular cavity. The performance of the method was quantified based on the L 2 norm error and the number of iterations needed to convergence the energy equation at each time step. For all three validation cases, a mesh convergence rate of approximately O(h) was obtained, which is below the expected accuracy of the numerical method. Only for the gallium melting case, the use of a higher-order method proved to be beneficial. The results from the present numerical campaign demonstrate the promise of the discontinuous Galerkin finite element method for modeling certain solid–liquid phase change problems where large gradients in the flow field are present or the phase change is highly localized, however, further enhancement of the method is needed to fully benefit from the use of a higher-order numerical method when solving solid–liquid phase change problems.

A fourth-order fractional Adams-type implicit–explicit method for nonlinear fractional ordinary differential equations with weakly singular solutions

Based on piecewise cubic interpolation polynomials, we develop a high-order numerical method with nonuniform meshes for nonlinear fractional ordinary differential equations (FODEs) with weakly singular solutions. This method is a fractional variant of the classical Adams-type implicit–explicit method widely used for integer order ordinary differential equations. We rigorously prove that the method is unconditionally convergent under the local Lipschitz condition of the nonlinear function, and when the mesh parameter is properly selected, it can achieve optimal fourth-order convergence for weakly singular solutions. We also prove the stability of the method, and discuss the applicability of the method to multi-term nonlinear FODEs and systems of multi-order nonlinear FODEs with weakly singular solutions. Numerical results are given to confirm the theoretical convergence results.

Summation-by-Parts Operators for General Function Spaces

Summation-by-parts (SBP) operators are popular building blocks for systematically developing stable and high-order accurate numerical methods for time-dependent differential equations. The main idea behind existing SBP operators is that the solution is assumed to be well approximated by polynomials up to a certain degree, and the SBP operator should therefore be exact for them. However, polynomials might not provide the best approximation for some problems, and other approximation spaces may be more appropriate. In this paper, a theory for SBP operators based on general function spaces is developed. We demonstrate that most of the established results for polynomial-based SBP operators carry over to this general class of SBP operators. Our findings imply that the concept of SBP operators can be applied to a significantly larger class of methods than is currently known. We exemplify the general theory by considering trigonometric, exponential, and radial basis functions.

Operator-Compensation Schemes Combining with Implicit Integration Factor Method for the Nonlinear Dirac Equation

A high-order accuracy numerical method for the (1+1)-dimensional nonlinear Dirac (NLD) equation is given in this work. For the spatial discretization, high-order operator-compensation technology is adopted, then semi-discrete scheme is obtained. Energy conservation and charge conservation are shown for the semi- discrete scheme. For the temporal discretization, implicit integration factor ( IIF) method is utilized to deal with the ordinary differential equations that are obtained from the semi-discrete scheme. The accuracy of the high-order numerical method is verified by numerical experiments, and the interaction dynamics of NLD solitary waves are investigated.

A hybrid finite element–finite volume method for conservation laws

We propose an arbitrarily high-order accurate numerical method for conservation laws that is based on a continuous approximation of the solution. The degrees of freedom are point values at cell interfaces and moments of the solution inside the cell. To lowest (3rd) order this method reduces to the Active Flux method. The update of the moments is achieved immediately by integrating the conservation law over the cell, integrating by parts and employing the continuity across cell interfaces. We propose two ways how the point values can be updated in time: either by first deriving a semi-discrete method that uses a finite-difference-type formula to approximate the spatial derivative, and integrating this method e.g. with a Runge-Kutta scheme, or by using a characteristics-based update, which is inspired by the original (fully discrete) Active Flux method. We analyze stability and accuracy of the resulting methods.

Experimental and numerical study of flame acceleration and DDT in a channel with continuous obstacles

Flame acceleration and deflagration-to-detonation transition (DDT) in obstructed channels is an important subject of research for propulsion and explosion safety. Experiment and numerical simulation of DDT in a stoichiometric hydrogen–oxygen mixture in a channel equipped with continuous triangular obstacles were conducted in this work. In the experiment, high-speed schlieren photography was used to record the evolution of reaction front and strong pressure waves. A pressure transducer was used to record the pressure build-up. In the numerical simulation, a high-order numerical method was used to solve the fully compressible reactive Navier–Stokes equations coupled with a calibrated chemical-diffusive model. The calculations are in good agreement with experimental observations. The result shows that the triangular obstacles can significantly promote flame acceleration and provide conditions for the occurrence of DDT. In the early stages of flame acceleration, the main cause for flame roll-up and distortion is the effect of vortices generated in the gaps between neighbouring triangular obstacles. The scales and velocities of vortices are determined by the positive feedback process between combustion-generated flow and flame propagation. The continuous triangular obstacles create an intricate flow field and increase the complexity of shock reflections. This complicated flow leads to local detonation initiation through different mechanisms, i.e. flame-flame collisions and flame-shock interactions. Successive local detonation ignitions and failures are produced in the obstacle gaps due to the continuous layout of the triangular obstacles. It was found that successive local detonation ignitions are critical for the eventual success of DDT formation because the shock waves generated by them continually strengthen the leading shock. The detonation failure or survival due to diffraction depends on the height of the narrow space (h*) between the bulk flame and obstacle vertex, and can be quantitatively characterised by the ratio of the space height to detonation cell size (), h*/.