High-order Time Integration Schemes Research Articles

Modeling gas flow in low-permeability formations: An efficient combination of mixed finite elements and high order time integration schemes

Numerical simulation of gas flow in low permeability formations is a challenging task due to the high nonlinearity induced by (i) the compressibility of the gas, (ii) the Klinkenberg slippage effect and (iii) the Langmuir adsorption of the gas on the pore surface. Because of these nonlinearities, modeling gas flow in low permeability formations requires a great deal of computational effort. In this work, we develop an efficient numerical model using advanced spatial and temporal discretization methods for a simultaneous solution of the coupled equations of gas flow, cubic Peng-Robinson equation of state, slippage effect and Langmuir adsorption. The spatial discretization is performed with the lumped hybrid formulation of the mixed finite element method which is well adapted for fluid flow in heterogeneous porous media. The time integration is performed with high-order methods via the method of lines (MOL) which allows large time steps and efficient solution of the nonlinear system of equations. Numerical experiments, performed for gas extraction in a heterogeneous domain, point out the high efficiency of the new model since it can be until 10 times more efficient than the classical first-order time discretization method. Results of global sensitivity analysis show that the intrinsic gas-phase permeability and the Klinkenberg factor are the most influential parameters controlling the pumping rate in the case of gas extraction in a homogeneous domain.

High-order composite implicit time integration schemes based on rational approximations for elastodynamics

A new approach for developing implicit composite time integration schemes is established starting with rational approximations of the matrix exponential in the solution of the equations of motion. The rational approximations are designed to have the same effective stiffness matrix in all sub-steps. An efficient algorithm is devised so that the implicit equation for a sub-step is in the same form as that of the trapezoidal method. The proposed M-schemes are Mth order accurate using M sub-steps. The amount of numerical dissipation is controlled by a user-specified parameter ρ∞. The (M+1)-schemes that are (M+1)th order accurate using M sub-steps can also be constructed, but the amount of numerical dissipation is built in the schemes and not adjustable. The order of accuracy of all the schemes is not affected by external forces and physical damping. Numerical examples demonstrate that the proposed high-order composite schemes are effective for introducing numerical dissipation and allow the use of large time step sizes in wave propagation problems. The source code written in MATLAB is available for download at: https://github.com/ChongminSong/CompsiteTimeIntegration.

Global high-order numerical schemes for the time evolution of the general relativistic radiation magneto-hydrodynamics equations

Modeling correctly the transport of neutrinos is crucial in some astrophysical scenarios such as core-collapse supernovae and binary neutron star mergers. In this paper, we focus on the truncated-moment formalism, considering only the first two moments (M1 scheme) within the grey approximation, which reduces Boltzmann seven-dimensional equation to a system of equations closely resembling the hydrodynamic ones. Solving the M1 scheme is still mathematically challenging, since it is necessary to model the radiation-matter interaction in regimes where the evolution equations become stiff and behave as an advection-diffusion problem. Here, we present different global, high-order time integration schemes based on Implicit-Explicit Runge-Kutta methods designed to overcome the time-step restriction caused by such behavior while allowing us to use the explicit Runge-Kutta commonly employed for the magneto-hydrodynamics and Einstein equations. Finally, we analyze their performance in several numerical tests.

High-order implicit time integration scheme with controllable numerical dissipation based on mixed-order Padé expansions

A single-step high-order implicit time integration scheme with controllable numerical dissipation at high frequencies is presented for the transient analysis of structural dynamic problems. The amount of numerical dissipation is controlled by a user-specified value of the spectral radius ρ∞ in the high frequency limit. Using this user-specified parameter as a weight factor, a Padé expansion of the matrix exponential solution of the equation of motion is constructed by mixing the diagonal and sub-diagonal expansions. An efficient time-stepping scheme is designed where systems of equations, similar in complexity to the standard Newmark method, are solved recursively. It is shown that the proposed high-order scheme achieves high-frequency dissipation, while minimizing low-frequency dissipation and period errors. The effectiveness of the provided dissipation control and the efficiency of the scheme are demonstrated by numerical examples. A simple guideline for the choice of the controlling parameter and time step size is provided. The source codes written in MATLAB and FORTRAN are available for download at: https://github.com/ChongminSong/HighOrderTimeIntegration.

Energy stability of exponential time differencing schemes for the nonlocal Cahn‐Hilliard equation

AbstractThe nonlocal Cahn‐Hilliard equation has attracted much attention these years. Despite the advantage of describing more practical phenomena for modeling phase transitions of microstructures in materials, the nonlocal operator in the equation brings a lot of extra computational costs compared with the local Cahn‐Hilliard equation. Thus high order time integration schemes are needed in numerical simulations. In this paper, we propose two classes of exponential time differencing (ETD) schemes for solving the nonlocal Cahn‐Hilliard equation. We first use the Fourier collocation method to discretize the spatial domain, and then the ETD‐based multistep and Runge‐Kutta schemes are adopted for the time integration. In particular, some specific multistep and Runge‐Kutta schemes up to fourth order are constructed. We rigorously establish the energy stabilities of the multistep schemes up to fourth order and the second order Runge‐Kutta scheme, which show that the first order ETD and the second order Runge‐Kutta schemes unconditionally decrease the original energy. We also theoretically prove the mass conservations of the proposed schemes. Several numerical experiments in two and three dimensions are carried out to test the temporal convergence rates of the schemes and to verify their mass conservations and energy stabilities. The long time simulations of coarsening dynamics are also performed to verify the power law for the energy decay.

High-Order Accurate Solutions of the Point Kinetics Equations with the Spectral Deferred Correction Method

Solving initial value problems with high-order methods receives considerable attention in many fields because these methods can potentially improve the accuracy of the simulation results with lower computational cost than low-order methods. Most methods, however, are either complicated to implement or unstable when the order of accuracy is high. The spectral deferred correction (SDC) method is a stable, robust, and efficient high-order time-integration scheme capable of an arbitrary order of accuracy. In this paper, we apply the SDC method to solve the initial value problem of the point kinetics equations (PKEs). For our implementation, we show that SDC is -stable for orders up to eight and the order of accuracy is verified for PKE problems with a range of different reactivities. A fifth-order SDC method was then implemented to solve the exact PKE in the transient multilevel method of MPACT. The error from solutions of the exact PKE with SDC is shown to be negligible. The investigations made here can provide the foundation for future investigations simulating the neutron transport problem using the high-order methods for both spatial discretization and time integration.

An efficient SBFEM-based approach for transient exterior vibro-acoustic analysis of power-law functionally graded shells

An efficient approach for transient exterior vibro-acoustic analysis of power-law functionally graded (FG) shells is developed using the scaled boundary finite element method. In the structural formulation, a shell element is treated as a three-dimensional continuum and its middle surface is represented with a quadrilateral spectral element. Along the thickness, the middle surface is scaled and the displacements are approximated using a quadratic Lagrange interpolation. The assumed natural strain method is applied to treat numerical locking and the integral along the thickness is performed analytically. In the acoustic formulation, velocity potential acts as the basic unknown. An infinite fluid is truncated by a spherical surface. The exterior field is simulated through the improved doubly asymptotic open boundary while the interior region is split into subdomains with their impedance evaluated by the improved continued-fraction method. The structure and fluid are discretized independently. The collocation procedure is utilized to perform the data transfer across the nonconforming interface. A high-order implicit time integration scheme is applied to solve the coupled system of equations. Numerical examples demonstrate that the proposed approach is stable, accurate and highly efficient. The effects of the geometric and material parameters on the vibro-acoustic behaviors of FG shells are systematically studied.

High-order methods for thermal behavior of plate structures with internal heat sources

The paper proposes solving transient heat transfer in plates using high-order isogeometric analysis and high-order time integration schemes. The problem is often faced in fire-structure interaction where the heat transfer is coupled with the stress analysis. A major advantage for the proposed approach comes from high order continuity between elements. Thus, the need for rotational degrees of freedom is eliminated when analyzing the stresses in plates and the same mesh can also be used to recover the heat transfer patterns. To achieve the full potential of such high-order finite elements in space, we also use high-order time integration schemes so that the numerical solution can be significantly more accurate than the standard approaches. Furthermore, the isogeometric analysis is also used to represent the exact geometry thanks to the basis functions generated from Non-Uniform Rational B-Splines. Several test examples for the transient diffusion problem are presented. Compared to standard methods and for a prescribed accuracy the proposed approach requires significantly fewer degrees of freedom and a corresponding improvement in the computational efficiency.

An SPH-based fully-Lagrangian meshfree implicit FSI solver with high-order discretization terms

In this paper, a novel SPH (Smoothed Particle Hydrodynamics)-based implicit structure model as well as its refined coupled fully-Lagrangian meshfree FSI (Fluid-Structure Interaction) solver are presented. The implicit structure model is founded through consideration of conservation law of linear momentum. The implicit time integration is achieved along with the assumption of linear elastic material. For development of refined implicit structure model, three novel schemes are considered; (1) High-order implicit Time integration (HT) scheme is considered for enhanced implicit update of particles’ velocities and positions; (2) Predictor-Corrector (PC) method is incorporated for improved estimation of rotation matrix; (3) High-order Discretization (HD) operator is employed for consistent calculation of deformation gradient and acceleration by elastic force. The proposed SPH-based implicit structure model is coupled with a refined ISPH fluid model so as to develop an enhanced fully-Lagrangian meshfree FSI solver, where consistent time coupling between incompressible fluid flow and stiff structures can be achieved by using the same time step sizes in both phases. The proposed implicit structure model as well as its corresponding FSI solver are verified in terms of robustness and accuracy by reproducing a set of well-known benchmark tests.

High-order implicit time integration scheme based on Padé expansions

A single-step high-order implicit time integration scheme for the solution of transient as well as wave propagation problems is presented. It is constructed from the Padé expansion of the matrix exponential solution of a system of first-order ordinary differential equations formulated in the state–space. By exploiting the techniques of polynomial factorization and partial fractions of rational functions, and by decoupling the solution for the displacement and velocity vectors a computationally efficient scheme is developed. An important feature of the novel algorithm is that no direct inversion of the mass matrix is required. Based on the diagonal Padé expansion of order M a time-stepping scheme of order 2M is constructed. When M is odd, one system of equation with a real matrix and (M−1)/2 systems of equations with complex matrices are solved recursively. When M is even, there are M/2 systems of equations with complex matrices. These systems are sparse and comparable in complexity to the standard Newmark method, i.e., the effective system matrix is a linear combination of the static stiffness, damping, and mass matrices. The derived second-order scheme is analytically equivalent to Newmark constant average acceleration method. The proposed time integrator has been implemented in MATLAB and FORTRAN using direct linear equation solvers. In this article, numerical examples featuring nearly one million degrees of freedom are presented to highlight the exceptional accuracy and efficiency of the novel time-stepping scheme in comparison with established second-order methods. The MATLAB- and FORTRAN-implementations are available from the authors upon request.

Effect of time integration scheme in the numerical approximation of thermally coupled problems: From first to third order

The advantages of using high-order time integration schemes for thermally coupled flows are assessed numerically. First-, second-, and third-order backward difference schemes are evaluated. The problem is solved in a decoupled manner using a nested iterative algorithm for the Navier–Stokes and energy equations to eliminate decoupling errors. For the space discretization, a stabilized finite element formulation of the variational multiscale type is applied to enable the use of equal order interpolation between the problem unknowns and ensure stable solutions for convection-dominated cases. The integration schemes are compared by solving the flow over a confined square including mixed heat convection in two and three dimensions. Improved numerical approximation of dynamic solutions using high-order schemes is demonstrated in the Richardson number range of 0≤|Ri|≤10 up to a Reynolds number of Re=225.

Maximum bound principle preserving integrating factor Runge–Kutta methods for semilinear parabolic equations

A large class of semilinear parabolic equations satisfy the maximum bound principle (MBP) in the sense that the time-dependent solution preserves for any time a uniform pointwise bound imposed by its initial and boundary conditions. The MBP plays a crucial role in understanding the physical meaning and the wellposedness of the mathematical model. Investigation on numerical algorithms with preservation of the MBP has attracted increasingly attentions in recent years, especially for the temporal discretizations, since the violation of MBP may lead to nonphysical solutions or even blow-ups of the algorithms. In this paper, we study high-order MBP-preserving time integration schemes by means of the integrating factor Runge–Kutta (IFRK) method. Beginning with the space-discrete system of semilinear parabolic equations, we present the IFRK method in general form and derive the sufficient conditions for the method to preserve the MBP. In particular, we show that the classic four-stage, fourth-order IFRK scheme is MBP preserving for some typical semilinear systems although not strong stability preserving, which can be instantly applied to the Allen–Cahn type of equations. To our best knowledge, this is the first time to present a fourth-order linear numerical method preserving the MBP. In addition, convergence of these numerical schemes is proved theoretically and verified numerically, as well as their efficiency by simulations of 2D and 3D long-time evolutional behaviors. Numerical experiments are also carried out for a model which is not a typical gradient flow as the Allen–Cahn type of equations.

Adaptive High-Order Fluid-Structure Interaction Simulations with Reduced Mesh-Motion Errors

This paper demonstrates an adaptive approach for solving fluid–structure interaction problems using high-fidelity numerical methods along with a detailed analysis of mesh-motion errors. A high-order partitioned approach is applied to couple the fluid and the structural subsystems, where the fluid subsystem is discretized using a discontinuous Galerkin finite-element method, while the structures solver uses a continuous Galerkin discretization. An explicit mapping is used as the primary mesh deformation algorithm. High-order time-integration schemes are used by the coupled solver to march forward in time, and the space-time mesh of the fluid subsystem is adapted using output-based methods. The error estimates for the unsteady outputs are evaluated by calculating the uncoupled, unsteady adjoint of the fluid subsystem. Adaptive meshing is used to demonstrate the importance of mesh-motion errors on output convergence, and a comprehensive analysis is conducted to control such errors arising from the mesh deformation algorithm. The benefits of adaptive meshing are demonstrated on a cantilevered Euler–Bernoulli beam placed in a low-Reynolds number flow and on a two-dimensional pitching–plunging airfoil in a high-Reynolds number flow.

Parallel Implicit-Explicit General Linear Methods

High-order discretizations of partial differential equations (PDEs) necessitate high-order time integration schemes capable of handling both stiff and nonstiff operators in an efficient manner. Implicit-explicit (IMEX) integration based on general linear methods (GLMs) offers an attractive solution due to their high stage and method order, as well as excellent stability properties. The IMEX characteristic allows stiff terms to be treated implicitly and nonstiff terms to be efficiently integrated explicitly. This work develops two systematic approaches for the development of IMEX GLMs of arbitrary order with stages that can be solved in parallel. The first approach is based on diagonally implicit multistage integration methods (DIMSIMs) of types 3 and 4. The second is a parallel generalization of IMEX Euler and has the interesting feature that the linear stability is independent of the order of accuracy. Numerical experiments confirm the theoretical rates of convergence and reveal that the new schemes are more efficient than serial IMEX GLMs and IMEX Runge-Kutta methods.

Efficient higher-order cycle jump integration of a continuum fatigue damage model

Simulating high-cycle fatigue with continuum models offers the possibility to model stress-redistributions, consider 3D-stress states and simplifies extensions to multi-physics problems. The computational cost of conventional cycle-by-cycle time integrations is reduced by reformulating the fatigue problem as an ordinary differential equation for the material state and solving it with high-order adaptive time integration schemes. The computational cost of calculating the change of the material state in one cycle is further reduced by a high-order fatigue-specific time integration. The approach is exemplarily demonstrated for a fatigue extension of the implicit gradient-enhanced damage model in 3D and compared to experimental Wöhler lines.

Robustness and efficiency of an implicit time-adaptive discontinuous Galerkin solver for unsteady flows

High-fidelity fluid dynamics simulations of unsteady flows are nowadays of great interest for many industrial fields. This class of simulations, as they are characterized by a wide range of temporal scales, requires robust, accurate and efficient long time integration strategies. These features can be achieved by an appropriate coupling of high-order time integration schemes and time-step adaptation algorithms. The adaptation algorithms are typically based on a local error estimator, which exploits the local truncation error of the time integration scheme and of its lower order embedded scheme.In literature few information are available to assess the benefits in terms of robustness, accuracy, and efficiency provided by the coupling between temporal schemes and adaptation strategies for unsteady CFD simulations. The aim of this work is to reduce this gap, presenting a numerical investigation of the performance for different adaptive time-step strategies, based on implicit Rosenbrock-type temporal schemes, in a high-order discontinuous Galerkin solver.The performance of the considered time integration strategies for the autonomous ODE system resulting from the DG space discretization of the Navier–Stokes equations is assessed for several test cases of increasing stiffness and difficulty, identifying the best scheme and algorithm: (i) the 2D laminar flow around a circular cylinder and around a tandem of cylinders at ReD=100; (ii) the 2D viscous flow through a porous media, modelled as an array of cylinders, at ReD=2100 and ReD=10,000; (iii) the 3D turbulent flow through a 4-wheels rudimentary landing gear (RLG) at ReD=1×106.

Fully implicit time integration in truly incompressible SPH

In the last years, the Smoothed Particle Hydrodynamics (SPH) method was introduced in new fields of engineering applications where highly viscous flow like polymer flow or very small geometries, like in micro-flow, are involved. In most of these applications it was not possible to use realistic fluid parameters, larger simulation domains or higher resolution because of restriction of time step due to the viscous time step criterion of the explicit time integration scheme. The computational effort was too high even for highly scalable SPH codes. In this article, we present a first-order implicit time stepping scheme for truly incompressible SPH (ISPH) to eliminate the viscous time step criterion. We propose a consistent time stepping approach where both velocity and particle position are solved implicitly and compare the results to traditional semi-implicit ISPH and implicit schemes where only velocity is integrated implicitly. We study Poiseuille flow, Taylor–Green flow and Rayleigh–Taylor instability and compare results of the implicit schemes to ISPH. Energy conservation is investigated using Taylor–Green vortex to highlight differences between the approaches. We find that the fully implicit time integration scheme is numerically stable. Since no error estimation of accuracy is used, the larger time step size leads to deviations of the trajectories. In future work, higher-order time integration schemes as well as error estimators should be investigated.

Implicit-solvent coarse-grained modeling for polymer solutionsviaMori-Zwanzig formalism

We present a bottom-up coarse-graining (CG) method to establish implicit-solvent CG modeling for polymers in solution, which conserves the dynamic properties of the reference microscopic system. In particular, tens to hundreds of bonded polymer atoms (or Lennard-Jones beads) are coarse-grained as one CG particle, and the solvent degrees of freedom are eliminated. The dynamics of the CG system is governed by the generalized Langevin equation (GLE) derived via the Mori-Zwanzig formalism, by which the CG variables can be directly and rigorously linked to the microscopic dynamics generated by molecular dynamics (MD) simulations. The solvent-mediated dynamics of polymers is modeled by the non-Markovian stochastic dynamics in GLE, where the memory kernel can be computed from the MD trajectories. To circumvent the difficulty in direct evaluation of the memory term and generation of colored noise, we exploit the equivalence between the non-Markovian dynamics and Markovian dynamics in an extended space. To this end, the CG system is supplemented with auxiliary variables that are coupled linearly to the momentum and among themselves, subject to uncorrelated Gaussian white noise. A high-order time-integration scheme is used to solve the extended dynamics to further accelerate the CG simulations. To assess, validate, and demonstrate the established implicit-solvent CG modeling, we have applied it to study four different types of polymers in solution. The dynamic properties of polymers characterized by the velocity autocorrelation function, diffusion coefficient, and mean square displacement as functions of time are evaluated in both CG and MD simulations. Results show that the extended dynamics with auxiliary variables can construct arbitrarily high-order CG models to reproduce dynamic properties of the reference microscopic system and to characterize long-time dynamics of polymers in solution.

Hybrid methods for radiation transport using diagonally implicit Runge–Kutta and space–time discontinuous Galerkin time integration

In this work, we describe extensions of a hybrid method for time-dependent linear, kinetic radiation transport problems to high-order time integration schemes of the diagonally-implicit Runge–Kutta (DIRK) and space–time discontinuous Galerkin (STDG) types. The hybrid methods are constructed by splitting the radiation flux into collided and uncollided components to which low- and high-resolution discrete ordinates approximations are applied. The efficiency of hybrid methods constructed using DIRK, STDG, integral deferred correction (IDC), and implicit Euler schemes is compared using a test problem in one-dimensional slab geometry containing material discontinuities. It is observed that (i) higher-order methods are more efficient than the implicit Euler method, often by an order of magnitude or more; (ii) third-order methods yield solutions of a given error in roughly half the time of the second-order method of the same type; and (iii) for a given order of accuracy it is found that the most efficient class of time integration scheme is STDG, followed by DIRK, with IDC the least efficient, for the test problem considered. Two test problems in two-dimensional xy-geometry are used to compare the computational efficiency of hybrid and standard discrete ordinates methods constructed with DIRK and STDG integrators. We observe that replacing a standard discrete ordinates method using an angular quadrature of order N with a hybrid discrete ordinates method using angular quadratures of order 2N and N/2 for the uncollided and collided fluxes, respectively, usually reduces overall solution time by a factor of 2 or more while simultaneously reducing the resulting solution error by a factor of 2 or more for the test problems considered.

Using finite element codes as a numerical platform to run molecular dynamics simulations

A mathematically rigorous methodology for embedding the governing equations of molecular dynamics in the formalism of the finite element method is presented. Only one generalized finite element type is needed to cover all different types of existing interatomic potentials. The finite element type is simply specified by two parameters characterizing the type of the interatomic potential to be considered. Built on this formulation a partitioned Runge–Kutta method—summarizing a wide range of explicit and implicit, single- and multi-stage, lower and higher order time integration schemes—is embedded in a unified manner. The required finite element residual vector and the related Jacobian matrix are stated explicitly. The related FE-mesh coincides with the neighborhood lists used in standard molecular dynamics enabling the use of common tools. The range, versatility and performance of the proposed finite element formulation have been demonstrated by means of several numerical examples.