High-order Accurate Schemes Research Articles

High-order accurate implicit-explicit time-stepping schemes for wave equations on overset grids

New implicit and implicit-explicit time-stepping methods for the wave equation in second-order form are described with application to two and three-dimensional problems discretized on overset grids. The implicit schemes are single step, three levels in time, and based on the modified equation approach. Second and fourth-order accurate schemes are developed and they incorporate upwind dissipation for stability on overset grids. The fully implicit schemes are useful for certain applications such as the WaveHoltz algorithm for solving Helmholtz problems where very large time-steps are desired. Some wave propagation problems are geometrically stiff due to localized regions of small grid cells, such as grids needed to resolve fine geometric features, and for these situations the implicit time-stepping scheme is combined with an explicit scheme: the implicit scheme is used for component grids containing small cells while the explicit scheme is used on the other grids such as background Cartesian grids. The resulting partitioned implicit-explicit scheme can be many times faster than using an explicit scheme everywhere. The accuracy and stability of the schemes are studied through analysis and numerical computations.

High-order accurate well-balanced energy stable finite difference schemes for multi-layer shallow water equations on fixed and adaptive moving meshes

This paper develops high-order accurate well-balanced (WB) energy stable (ES) finite difference schemes for multi-layer (the number of layers M⩾2) shallow water equations (SWEs) with non-flat bottom topography on both fixed and adaptive moving meshes, extending our previous work on the single-layer shallow water magnetohydrodynamics [25] and single-layer SWEs on adaptive moving meshes [58]. To obtain an energy inequality, the convexity of an energy function for an arbitrary M is proved by finding recurrence relations of the leading principal minors or the quadratic forms of the Hessian matrix of the energy function with respect to the conservative variables, which is more involved than the single-layer case due to the coupling between the layers in the energy function. An important ingredient in developing high-order semi-discrete ES schemes is the construction of a two-point energy conservative (EC) numerical flux. In pursuit of the WB property, a sufficient condition for such EC fluxes is given with compatible discretizations of the source terms similar to the single-layer case. It can be decoupled into M identities individually for each layer, making it convenient to construct a two-point EC flux for the multi-layer system. To suppress possible oscillations near discontinuities, WENO-based dissipation terms are added to the high-order WB EC fluxes, which gives semi-discrete high-order WB ES schemes. Fully-discrete schemes are obtained by employing high-order explicit strong stability preserving Runge-Kutta methods and proved to preserve the lake at rest. The schemes are further extended to moving meshes based on a modified energy function for a reformulated system, relying on the techniques proposed in [58]. Numerical experiments are conducted for some two- and three-layer cases to validate the high-order accuracy, WB and ES properties, and high efficiency of the schemes, with a suitable amount of dissipation chosen by estimating the maximal wave speed due to the lack of an analytical expression for the eigenstructure of the multi-layer system.

An efficient preconditioner for linear systems arising from high-order accurate schemes of time fractional diffusion equations

An efficient preconditioner for linear systems arising from high-order accurate schemes of time fractional diffusion equations

Stabilizing Discontinuous Galerkin Methods Using Dafermos’ Entropy Rate Criterion: II—Systems of Conservation Laws and Entropy Inequality Predictors

A novel approach for the stabilization of the Discontinuous Galerkin method based on the Dafermos entropy rate crition is presented. First, estimates for the maximal possible entropy dissipation rate of a weak solution are derived. Second, families of conservative Hilbert–Schmidt operators are identified to dissipate entropy. Steering these operators using the bounds on the entropy dissipation results in high-order accurate shock-capturing DG schemes for the one-dimensional Euler equations, satisfying the entropy rate criterion and an entropy inequality. Other testcases include the one-dimensional Buckley–Leverett equation.

Convergence towards High-Speed Steady States Using High-Order Accurate Shock-Capturing Schemes

Creating time-marching unsteady governing equations for a steady state in high-speed flows is not a trivial task. Residue convergence in time cannot be achieved when using most low- and high-order spatial discretization schemes. Recently, high-order, weighted, essentially non-oscillatory schemes have been specially designed for steady-state simulations. They have been shown to be capable of achieving machine precision residues when simulating the Euler equations under canonical coordinates. In the present work, we review these schemes and show that they can also achieve machine residues when simulating the Navier–Stokes equations under generalized coordinates. This is carried out by considering three supersonic flows of perfect fluids, namely the flow upstream a cylinder, the flow over a blunt wedge, and the flow over a compression ramp.

High-order discretization–based self-adaptive turbulence eddy simulation for supersonic base flow with PHengLEI software

PurposeThe purpose of the present study is to develop a new numerical framework that can predict the supersonic base flow more accurately, including the development of axisymmetrically separated shear layer and recompression shock. To this end, two aspects are improved and combined, i.e. a newly self-adaptive turbulence eddy simulation (SATES) turbulence modeling method and a high-order discretization numerical scheme. Furthermore, the performance of the new numerical framework within a general-purpose PHengLEI software is assessed in detail.Design/methodology/approachSatisfactory prediction of the supersonic separated shear layer with unsteady wake flow is quite challenging. By using a unified turbulence model called SATES combining high-order accurate discretization numerical schemes, the present study first assesses the performance of newly developed SATES for supersonic axisymmetric separation flows. A high-order finite differencing-based compressible computational fluid dynamics (CFD) code called PHengLEI is developed and several different numerical schemes are used to investigate the effects on shock-turbulence interactions, which include the monotonic upstream-centered scheme for conservation laws (MUSCL), weighted compact nonlinear scheme (WCNS) and hybrid cell-edge and cell-node dissipative compact scheme (HDCS).FindingsCompared with the available experimental data and the numerical predictions, the results of SATES by using high-order accurate WCNS or HDCS schemes agree better with the experiments than the results by using the MUSCL scheme. The WCNS and HDCS can also significantly improve the prediction of flow physics in terms of the instability of the annular shear layer and the evolution of the turbulent wake.Research limitations/implicationsThe small deviations in the recirculation region can be found between the present numerical results and experimental data, which could be caused by the inaccurate incoming boundary layer condition and compressible effects. Therefore, a proper incoming boundary layer condition with turbulent fluctuations and compressibility effects need to be considered to further improve the accuracy of simulations.Practical implicationsThe present study evaluates a high-order discretization-based SATES turbulence model for supersonic separation flows, which is quite valuable for improving the calculation accuracy of aeronautics applications, especially in supersonic conditions.Originality/valueFor the first time, the newly developed SATES turbulence modeling method combining the high-order accurate WCNS or HDCS numerical schemes is implemented on the PHengLEI software and successfully applied for the simulations of supersonic separation flows, and satisfactory results are obtained. The unsteady evolutions of the supersonic annular shear layer are analyzed, and the hairpin vortex structures are found in the simulation.

High-order accurate variable time step compact schemes for pricing vanilla and exotic options

High-order accurate variable time step compact schemes for pricing vanilla and exotic options

Temporal High-Order Accurate Numerical Scheme for the Landau–Lifshitz–Gilbert Equation

In this paper, a family of temporal high-order accurate numerical schemes for the Landau–Lifshitz–Gilbert (LLG) equation is proposed. The proposed schemes are developed utilizing the Gauss–Legendre quadrature method, enabling them to achieve arbitrary high-order time discretization. Furthermore, the geometrical properties of the LLG equation, such as the preservation of constant magnetization magnitude and the Lyapunov structure, are investigated based on the proposed discrete schemes. It is demonstrated that the magnetization magnitude remains constant with an error of (2p+3) order in time when utilizing a (2p+2)th-order discrete scheme. Additionally, the preservation of the Lyapunov structure is achieved with a second-order error in the temporal step size. Numerical experiments and simulations effectively verify the performance of our proposed algorithm and validate our theoretical analysis.

Fully decoupled and high-order linearly implicit energy-preserving RK-GSAV methods for the coupled nonlinear wave equation

This paper is concerned with high-order accurate, linearly implicit and energy-preserving schemes for the coupled nonlinear wave equation. To this end, a novel auxiliary variable approach proposed in a recent paper (Ju et al., 2022) is applied to reformulate the governing equation in an equivalent system, which possesses a modified energy that consists of primary functional and quadratic functional. Subsequently, we utilize the extrapolation strategy/prediction–correction technique for treating the nonlinear terms of the equivalent system to obtain a linearized energy-preserving system. Then, a series of fully decoupled high-order linear energy-preserving schemes are constructed by using the symplectic Runge–Kutta (RK) methods. We show that the modified energy functional can be precisely conserved under certain circumstances for the coefficients of the symplectic RK methods, and the unique solvability of the proposed time-stepping scheme is analyzed. In numerical implementations, a class of linear algebraic equations with constant-coefficient matrices need to be solved at each time step and only less computational costs are required. Extensive numerical experiments of the coupled nonlinear wave equation with local/nonlocal diffusion operators are provided to illustrate the energy conservation law and the computational efficiency of the schemes in long-time simulations.

High-Order Accurate Entropy Stable Schemes for Relativistic Hydrodynamics with General Synge-Type Equation of State

High-Order Accurate Entropy Stable Schemes for Relativistic Hydrodynamics with General Synge-Type Equation of State

Computation of unsteady electromagnetic scattering about 3D complex bodies in free space with high-order difference potentials

We extend the previously developed high-order accurate method for acoustic scattering to electromagnetic scattering, i.e., from the scalar setting to a vector setting. First, the governing Maxwell's equations are reduced from their original first-order form to a system of second-order wave equations for the individual Cartesian components of electromagnetic field. In free space, these wave equations are uncoupled. Yet at the boundary of the scatterer, the variables that they govern (i.e., Cartesian field components) remain fully coupled via the boundary conditions that account for the specific scattering mechanism. Next, the wave equations are equivalently replaced with Calderon's boundary equations with projections obtained using the method of difference potentials and a compact high-order accurate scheme. The Calderon's boundary equations are combined with the boundary conditions and the overall system is solved by least squares. The resulting vector methodology (electromagnetic) inherits many useful properties of the scalar one (acoustic). In particular, it offers sub-linear computational complexity, does not require any special treatment of the artificial outer boundary, and has the capacity to solve multiple similar problems at a low individual cost per problem. We demonstrate the performance of the new method by computing the scattering of a given impinging wave about a double-cone hypersonic shape.

Study on bubble dynamics in sodium using three-dimensional MPS method

Sodium-cooled fast reactor (SFR) is one of the most promising Generation IV reactors. Design and safety evaluation of a SFR requires sufficient knowledge about gas bubble behavior in sodium pool. It is challenging to perform numerical study on bubble dynamics in sodium pool due to the discontinuity caused by the differences of physical properties. In this study, a three-dimensional (3D) Moving Particle Semi-implicit (MPS) method is presented featuring high-order accurate scheme and Riemann-based stable multiphase model. Numerical discontinuity caused by large density, large viscosity ratio and surface tension can be handled by the proposed method. To verify the proposed method, two cases from the bubble diagram of Grace (1973) is simulated and rising velocities and shapes of rising gas bubbles are compared, both are agreeing very well. Then co-axial coalescence of gas bubbles are calculated and compared with existing experimental and numerical results. Finally bubble rising behavior in sodium is analyzed using the developed 3D MPS method.

Numerical simulation of acoustic streaming in standing waves

This study presents a numerical investigation of acoustic streaming motion (of the Rayleigh type) in a compressible gas inside two-dimensional rectangular enclosures. To numerically study the effects of the sound field intensity on the formation process of streaming structures, we propose to discretize the fully compressible form of the two-dimensional Navier-Stokes equations using a high-order (formally greater or equal to fourth-order) accurate numerical scheme in both space and time. The proposed numerical solver utilizes high-order compact schemes along each spatial dimension combined with a filtering procedure when it is necessary. Acoustic standing waves are excited inside the enclosures and the resulting acoustic streaming patterns are investigated for low and high-intensity waves, in both linear and nonlinear regimes following closely the work of Daru et al. (2013). An extended investigation indicates that without the incorporation of an appropriate filter, the application of the high-order compact schemes in case of fast streaming results in spurious oscillations which inhibit their applicability. Following the recent relevant literature, the numerical simulations performed demonstrate the ability of the proposed numerical approach to reproduce efficiently and robustly the transitions from regular acoustic streaming to irregular streaming and the relevant phenomena confirming results that have been previously presented.

Self-Adaptive Turbulence Eddy Simulation with a high-order finite differencing method for high Reynolds number complex flows

Accurate prediction of high Reynolds number wall-bounded flows is a challenge for turbulence modeling and numerical discretization. Moreover, complex engineering applications increase the difficulties in implementing robust and low-dissipation numerical schemes. A newly developed unified turbulence approach, Self-Adaptive Turbulence Eddy Simulation (SATES), is assessed in the simulation of a multi-element 30P30N airfoil at a high Reynolds number of 1.71 × 106. To reduce the numerical error, high-order finite differencing-based weighted compact nonlinear schemes (WCNSs) are successfully applied. The effects of grid resolution, numerical schemes and turbulence model are investigated by grouping studies. The SATES shows low sensitivity to grid resolution and predicts satisfactory results even on a coarse grid containing approximately 7.2 million cells. Compared with second-order accurate central differencing, the high-order accurate WCNS has less numerical dissipation and significantly improves the resolution of turbulence structures and the accuracy of numerical results. Z-type nonlinear weighted averaging interpolation, designed for shock capturing property in the WCNS scheme, leads to extra numerical dissipation and affects the accuracy of the results. The numerical results of the SATES are more accurate than those of the Smagorinsky Large Eddy Simulation (LES) and the Improved Delayed Detached Eddy Simulation (IDDES). The SATES provides sufficient modeling of the boundary layer on the coarse grid and a fast ‘transition’ from the RANS to the LES mode for accurate prediction of the separated shear layer and flap boundary layer separation. These results confirm that the SATES turbulence modeling approach is efficient for high Reynolds number flows and that the high-order accurate WCNS schemes show sufficient robustness and accuracy for engineering applications.

Imposing slip conditions on curved boundaries for 3D incompressible flows with a very high-order accurate finite volume scheme on polygonal meshes

The conventional no-slip boundary condition does not always hold in several fluid flow applications and must be replaced with appropriate slip conditions according to the wall and fluid properties. However, not only slip boundary conditions are still a subject of discussion among fluid dynamicists, but also their numerical treatment is far from being well-developed, particularly in the context of very high-order accurate methods. The complexity of these conditions significantly increases when the boundary is not aligned with the chosen coordinate system and, even more challenging, when the fluid slips along a curved boundary. The present work proposes a simple and efficient numerical treatment of general slip boundary conditions on arbitrary curved boundaries for three-dimensional fluid flow problems governed by the incompressible Navier–Stokes equations. In that regard, two critical challenges arise: (i) achieving very high-order of convergence with arbitrary curved boundaries for the classical no-slip boundary conditions and (ii) extending the developed numerical techniques to impose general slip boundary conditions. The conventional treatment of curved boundaries relies on generating curved meshes to eliminate the geometrical mismatch between the physical and computational boundaries and achieve high-order of convergence. However, such an approach requires sophisticated meshing algorithms, cumbersome quadrature rules on curved elements, and complex non-linear transformations. In contrast, the reconstruction for off-site data (ROD) method handles arbitrary curved boundaries approximated with linear piecewise elements while employing polynomial reconstructions with specific linear constraints to fulfil the prescribed boundary conditions. For that purpose, the general slip boundary conditions are reformulated on a local orthonormal basis to allow a straightforward application of the ROD method with scalar boundary conditions. The Navier–Stokes equations are discretised with a staggered finite volume method, and the numerical fluxes are computed solely on the polygonal mesh elements. Several benchmark test cases of fluid flow problems in non-trivial three-dimensional curved domains are addressed and confirm that the proposed method effectively achieves up to the eighth-order of convergence.

On the Convergence of α Schemes with Source Terms for Scalar Convex Conservation Laws

In this study, we use an extension of Yang’s convergence criterion [N. Jiang, On the wavewise entropy inequality for high-resolution schemes with source terms II: the fully discrete case] to show the entropy convergence of a class of fully discrete α schemes, now with source terms, for non-homogeneous scalar convex conservation laws in the one-dimensional case. The homogeneous counterparts (HCPs) of these schemes were constructed by S. Osher and S. Chakravarthy in the mid-1980s [A New Class of High Accuracy TVD Schemes for Hyperbolic Conservation Laws (1985), Very High Order Accurate TVD Schemes (1986)], and the entropy convergence of these methods, when m=2, was settled by the author [N. Jiang, The Convergence of α Schemes for Conservation Laws II: Fully-Discrete]. For semi-discrete α schemes, with or without source terms, the entropy convergence of these schemes was previously established (for m=2) by the author [N. Jiang, The Convergence of α Schemes for Conservation Laws I: Semi-Discrete Case].

A fast dynamic smooth adaptive meshing scheme with applications to compressible flow

We develop a fast-running smooth adaptive meshing (SAM) algorithm for dynamic curvilinear mesh generation, which is based on a fast solution strategy of the time-dependent Monge-Ampère (MA) equation, det⁡∇ψ(x,t)=G∘ψ(x,t). The novelty of our approach is a new so-called perturbation formulation of MA, which constructs the solution map ψ via composition of a sequence of near-identity deformations of a reference mesh. Then, we formulate a new version of the deformation method [21] that results in a simple, fast, and high-order accurate numerical scheme and a dynamic SAM algorithm that is of optimal complexity when applied to time-dependent mesh generation for solutions to hyperbolic systems such as the Euler equations of gas dynamics. We perform a series of challenging 2D and 3D mesh generation experiments for grids with large deformations, and demonstrate that SAM is able to produce smooth meshes comparable to state-of-the-art solvers [22,18], while running approximately 200 times faster. The SAM algorithm is then coupled to a simple Arbitrary Lagrangian Eulerian (ALE) scheme for 2D gas dynamics. Specifically, we implement the C-method [64,65] and develop a new ALE interface tracking algorithm for contact discontinuities. We perform numerical experiments for both the Noh implosion problem as well as a classical Rayleigh-Taylor instability problem. Results confirm that low-resolution simulations using our SAM-ALE algorithm compare favorably with high-resolution uniform mesh runs.

A new family of semi-implicit Finite Volume/Virtual Element methods for incompressible flows on unstructured meshes

We introduce a new family of high order accurate semi-implicit schemes for the solution of nonlinear time-dependent systems of partial differential equations (PDE) on unstructured polygonal meshes. The time discretization is based on a splitting between explicit and implicit terms that may arise either from the multi-scale nature of the governing equations, which involve both slow and fast scales, or in the context of projection methods, where the numerical solution is projected onto the physically meaningful solution manifold. We propose to use a high order finite volume (FV) scheme for the explicit terms, hence ensuring conservation property and robustness across shock waves, while the virtual element method (VEM) is employed to deal with the discretization of the implicit terms, which typically requires an elliptic problem to be solved. The numerical solution is then transferred via suitable L2 projection operators from the FV to the VEM solution space and vice-versa. High order time accuracy is then achieved using the semi-implicit IMEX Runge–Kutta schemes, and the novel schemes are proven to be asymptotic preserving (AP) and well-balanced (WB). As representative models, we choose the shallow water equations (SWE), thus handling multiple time scales characterized by a different Froude number, and the incompressible Navier–Stokes equations (INS), which are solved at the aid of a projection method to satisfy the solenoidal constraint of the velocity field. Furthermore, an implicit discretization for the viscous terms is also devised for the INS model, which is based on the VEM technique. Consequently, the CFL-type stability condition on the maximum admissible time step is based only on the fluid velocity and not on the celerity nor on the viscous eigenvalues. A large suite of test cases demonstrates the accuracy and the capabilities of the new family of schemes to solve relevant benchmarks in the field of incompressible fluids.

A New Entropy Stable Finite Difference Scheme for Hyperbolic Systems of Conservation Laws

The hyperbolic problem has a unique entropy solution, which maintains the entropy inequality. As such, people hope that the numerical results should maintain the discrete entropy inequalities accordingly. In view of this, people tend to construct entropy stable (ES) schemes. However, traditional numerical schemes cannot directly maintain discrete entropy inequalities. To address this, we here construct an ES finite difference scheme for the nonlinear hyperbolic systems of conservation laws. The proposed scheme can not only maintain the discrete entropy inequality, but also enjoy high-order accuracy. Firstly, we construct the second-order accurate semi-discrete entropy conservative (EC) schemes and ensure that the schemes meet the entropy identity when an entropy pair is given. Then, the second-order EC schemes are used as a building block to achieve the high-order accurate semi-discrete EC schemes. Thirdly, we add a dissipation term to the above schemes to obtain the high-order ES schemes. The term is based on the Weighted Essentially Non-Oscillatory (WENO) reconstruction. Finally, we integrate the scheme using the third-order Runge–Kutta (RK) approach in time. In the end, plentiful one- and two-dimensional examples are implemented to validate the capability of the scheme. In summary, the current scheme has sharp discontinuity transitions and keeps the genuine high-order accuracy for smooth solutions. Compared to the standard WENO schemes, the current scheme can achieve higher resolution.

High order accurate and convergent numerical scheme for the strongly anisotropic Cahn–Hilliard model

AbstractWe propose and analyze a second order accurate in time, energy stable numerical scheme for the strongly anisotropic Cahn–Hilliard system, in which a biharmonic regularization has to be introduced to make the equation well‐posed. A convexity analysis on the anisotropic interfacial energy is necessary to overcome an essential difficulty associated with its highly nonlinear and singular nature. The second order backward differentiation formula temporal approximation is applied, combined with Fourier pseudo‐spectral spatial discretization. The nonlinear surface energy part is updated by an explicit extrapolation formula. Meanwhile, the energy stability analysis is enforced by the fact that all the second order functional derivatives of the energy stay uniformly bounded by a global constant. A Douglas‐Dupont type regularization is added to stabilize the numerical scheme, and a careful estimate ensures a modified energy stability with a uniform constraint for the regularization parameter . In turn, the combination with an appropriate treatment for the nonlinear double well potential terms leads to a weakly nonlinear scheme. More importantly, such an energy stability is in terms of the interfacial energy with respect to the original phase variable, which enables us to derive an optimal rate convergence analysis.