Time Step Size Limitation Research Articles

An implicit–explicit relaxation extrapolated Runge–Kutta and energy-preserving finite element method for Klein–Gordon–Schrödinger equations

An implicit–explicit (IMEX) relaxation extrapolated Runge–Kutta (RERK) and energy-preserving finite element method is designed for Klein–Gordon–Schrödinger (KGS) equations with periodic boundary conditions. First, the RERK method is employed to the discretization in time, which can keep the unconditionally stable and energy conservation after selecting the appropriate relaxation parameters γn. An IMEX time-marching discretization is constructed so that only a system of linear equations needs to be solved at any time step. Next, the finite element method is utilized to the discretization in space. By the temporal–spatial splitting technique, optimal error estimates in the L2-norm and H1-norm are obtained with the convergent order O(τs+hk+1) and O(τs+hk), without any time stepsize restrictions. At last, some numerical examples are presented to demonstrate the theoretical results.

Uniformly accurate integrators for Klein–Gordon–Schrödinger systems from the classical to non-relativistic limit regime

In this paper we present a novel class of asymptotic consistent exponential-type integrators for Klein–Gordon–Schrödinger systems that capture all regimes from the slowly varying classical regime up to the highly oscillatory non-relativistic limit regime. We achieve convergence of order one and two that is uniform in c without any time step size restrictions. In particular, we establish an explicit relation between gain in negative powers of the potentially large parameter c in the error constant and loss in derivative.

A high-order generalized Finite Element Method for multiscale structural dynamics and wave propagation

This paper introduces a high-order multiscale Generalized/eXtended Finite Element Method (GFEM) tailored for the solution of structural dynamics and wave propagation problems exhibiting fine-scale and/or localized solution features such as singularities and discontinuities. The proposed method uses an explicit time-marching scheme together with a block-diagonal lumped mass matrix applicable to arbitrary patch approximation spaces and adopts shape functions computed numerically on-the-fly through the solution of local initial/boundary value problems. Numerical results in one-, two-, and three-dimensions for problems exhibiting singularities in the spatial gradient show that the proposed GFEM is able to capture relevant features of the response using structural-scale meshes that are much coarser than those required by Direct Generalized Finite Element Analyses (DGFEAs) of comparable accuracy. In addition, a detailed study of the critical time step size of the method is presented and compared against the one for DGFEA discretizations that provide similar levels of accuracy. It shows that the proposed method has considerably looser time step size restrictions than available DGFEAs.

Adjoint complement to the volume-of-fluid method for immiscible flows

The paper is concerned with an adjoint complement to the Volume-of-Fluid (VoF) method for immiscible two-phase flows, e.g. air and water, which is widely used in marine engineering due to its computational efficiency. The particular challenge of the primal and the corresponding adjoint VoF-approach refers to the sharp interface treatment featuring discontinuous physical properties. Both the continuous adjoint two-phase system (integration-by-parts) and the corresponding dual compressive convection schemes (summation-by-parts) are derived for two prominent compressive convection schemes, namely the High Resolution Interface Capturing Scheme (HRIC) and Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAM). The dual scheme rigorously mirrors the primal Normalized-Variable-Diagram (NVD) stencils. Attention is restricted to steady state applications. Thus both the primal and the dual procedures are performed in pseudo-time and the backward integration of the dual approach is performed around the (pseudo-temporal) converged primal field. Therefore, the adjoint system experiences the same time step size restrictions as the primal system, is independent of the primal time horizon and forms a robust as well as an a priori stable adjoint solution process.The paper analyses the primal and adjoint equations for an engineering model problem. An analytical solution to the model problem is initially presented, which displays that the adjoint part does not offer a unique, non-trivial solution. As a remedy, an additional diffusive concentration term is introduced to the adjoint concentration equation. The imposed heuristic modification violates the dual consistency but strongly regularizes the solution of the adjoint system. The modification can be justified by reference to phase-separating diffuse-interface models and inheres a free mobility-parameter. Numerical results obtained from the modified approach are benchmarked against the analytical solution for the model problem. Supplementary, the influence of the modification on the sensitivities obtained from simulations for the two-dimensional flow around a submerged hydrofoil at Froude and Reynolds numbers of practical interest are discussed for a range of mobility-parameters. The final application refers to a shape-optimization of a generic 3D underwater vehicle and underlines a negligible influence of the free mobility parameter, even for an objective functional that directly depends on the manipulated (dual) field quantity.

Implicit Material Point Method with Convected Particle Domain Interpolation for Consolidation and Dynamic Analysis of Saturated Porous Media with Massive Deformation

A coupling convected particle domain interpolation based implicit material point method (ICCPDI) is developed in this paper for the long-time consolidation and dynamic problems of fully saturated porous media involving massive deformation. In the method, the implicit formulas based on the u-p form governing equations are derived to overcome the time step size limitation within the conventional explicit method and in turn improve the computational efficiency particularly for consolidation problems. Then, the CPDI technique is adopted to eliminate the numerical noises occurring while material points crossing cell boundaries in the massive deformation of solid skeleton. Moreover, an extended material point penalty method is proposed to effectively treat the boundary constraints in the coupled system by avoiding introducing an auxiliary boundary layer. Numerical simulation results of several representative examples demonstrate the accuracy and efficiency of the proposed ICCPDI method to predict the consolidation and dynamic responses of saturated porous media with massive deformation.

Analysis of Lowest-Order Characteristics-Mixed FEMs for Incompressible Miscible Flow in Porous Media

The time discrete scheme of characteristics type is especially effective for convection-dominated diffusion problems. The scheme has been used in various engineering areas with different approximations in spatial direction. The lowest-order mixed method is the most popular one for miscible flow in porous media. The method is based on a linear Lagrange approximation to the concentration and the zero-order Raviart-Thomas approximation to the pressure/velocity. However, the optimal error estimate for the lowest-order characteristics-mixed FEM has not been presented although numerous effort has been made in last several decades. In all previous works, only first-order accuracy in spatial direction was proved under certain time-step and mesh size restrictions. The main purpose of this paper is to establish optimal error estimates, $i.e.$, the second-order in $L^2$-norm for the concentration and the first-order for the pressure/velocity, while the concentration is more important physical component for the underlying model. For this purpose, an elliptic quasi-projection is introduced in our analysis to clean up the pollution of the numerical velocity through the nonlinear dispersion-diffusion tensor and the concentration-dependent viscosity. Moreover, the numerical pressure/velocity of the second-order accuracy can be obtained by re-solving the (elliptic) pressure equation at a given time level with a higher-order approximation. Numerical results are presented to confirm our theoretical analysis.

Application of the generalized finite-element method to the acoustic wave simulation in exploration seismology

The generalized finite-element method (GFEM) has been applied frequently to solve harmonic wave equations, but its use in the simulation of transient wave propagation is still limited. We have applied GFEM for the simulation of the acoustic wave equation in models relevant to exploration seismology. We also perform an assessment of its accuracy and efficiency. The main advantage of GFEM is that it provides an enhanced solution accuracy compared to the standard finite-element method (FEM). This is attained by adding user-defined enrichment functions to standard FEM approximations. For the acoustic wave equation, we consider plane waves oriented in different directions as the enrichments, whose argument includes the largest wavenumber of the wavefield. We combine GFEM with an unconditionally stable time integration scheme with a constant time step. To assess the accuracy and efficiency of GFEM, we compare the GFEM simulation results against those obtained with the spectral-element method (SEM). We use SEM because it is the method of choice for wave propagation simulation due to its proven accuracy and efficiency. In the numerical examples, we first perform a convergence study in space and time, evaluating the accuracy of both methods against a semianalytical solution. Then, we consider simulations of relevant models in exploration seismology that include low-velocity features, an inclusion with a complex geometric boundary and topography. Results using these models indicate that GFEM presents comparable accuracy and efficiency to those based on SEM. For the given examples, the GFEM efficiency stems from the combined effect of local mesh refinement, nonconforming or unstructured, and the unconditionally stable time integration scheme with a constant time step. Moreover, these features provide great flexibility to the GFEM implementations, proving advantageous when using, for example, unstructured grids that impose severe time step size restrictions in SEM.

Time adaptive conservative finite volume method

Finite volume methods for time dependent problems typically employ time integration schemes with constant time step sizes. The latter, in order to ensure stability and time accurate solutions, have to be chosen small enough, such that the Courant-Friedrichs-Lewy (CFL) number stays below a critical value everywhere. In many cases, however, this time step size limitation leads to huge numbers of very small time steps, which renders simulations very expensive. Motivated by this drawback of conventional time stepping schemes, various sub-time stepping algorithms have been proposed. However, most of them are inherently asynchronous, require small local CFL numbers or are not strictly conservative. In this paper a new adaptive time integration scheme for finite volume methods, which is conservative, of high spatial and temporal order, robust and easy to implement is presented. It relies on local sub-time steps which are fractions of a global time step by powers of two, i.e., the grid cells proceed asynchronously in the order of their termination time, but since they all synchronize at the end of each global time step, it is possible to guarantee continuity of the mean fluxes and thus strict conservation at the global time step resolution. Numerical experiments with 1D and 2D test cases demonstrate that the adaptive conservative time integration (ACTI) scheme can achieve extreme speed-up factors over conventional time integration, while still maintaining high spatial and temporal accuracy.

An efficient second-order linear scheme for the phase field model of corrosive dissolution

In this paper we propose an efficient numerical scheme for solving the phase field model (PFM) of corrosive dissolution that is linear and second-order accurate in both time and space. The PFM of corrosion is based on the gradient flow of a free energy functional depending on a phase field variable and a single concentration variable. While classic backward differentiation formula (BDF) schemes have been used for time discretization in the literature, they require very small time step sizes owing to the strong numerical stiffness and nonlinearity of the parabolic partial differential equation (PDE) system defining the PFM. Based on the observation that the governing equation corresponding to the phase field variable is very stiff due to the reaction term, the key idea of this paper is to employ an exponential time integrator that is more effective for stiff dynamic PDEs. By combining the exponential integrator based Rosenbrock–Euler scheme with the classic Crank–Nicolson scheme for temporal integration of the spatially semi-discretized system, we develop a decoupled linear numerical scheme that alleviates the time step size restriction due to high stiffness. Several numerical examples are presented to demonstrate accuracy, efficiency and robustness of the proposed scheme in two-dimensions, and we find that a time step size of 10−3 second for meshes with the typical spatial resolution 1 μm is stable. Additionally, the proposed scheme is robust and does not suffer from any convergence issues often encountered by nonlinear Newton methods.

Randomized Hamiltonian Monte Carlo

Tuning the durations of the Hamiltonian flow in Hamiltonian Monte Carlo (also called Hybrid Monte Carlo) (HMC) involves a tradeoff between computational cost and sampling quality, which is typically challenging to resolve in a satisfactory way. In this article we present and analyze a randomized HMC method (RHMC), in which these durations are i.i.d. exponential random variables whose mean is a free parameter. We focus on the small time step size limit, where the algorithm is rejection-free and the computational cost is proportional to the mean duration. In this limit, we prove that RHMC is geometrically ergodic under the same conditions that imply geometric ergodicity of the solution to underdamped Langevin equations. Moreover, in the context of a multi-dimensional Gaussian distribution, we prove that the sampling efficiency of RHMC, unlike that of constant duration HMC, behaves in a regular way. This regularity is also verified numerically in non-Gaussian target distributions. Finally we suggest variants of RHMC for which the time step size is not required to be small.

Fast fluid–structure interaction simulations using a displacement-based finite element model equipped with an explicit streamline integration prediction

We propose here a displacement-based updated Lagrangian fluid model developed to facilitate a monolithic coupling with a wide range of structural elements described in terms of displacements. The novelty of the model consists in the use of the explicit streamline integration for predicting the end-of-step configuration of the fluid domain. It is shown that this prediction considerably alleviates the time step size restrictions faced by the former Lagrangian models due to the possibility of an element inversion within one time step. The method is validated and compared with conventional approaches using three numerical examples. Time step size and corresponding Courant numbers leading to optimal behavior in terms of computational efficiency are identified.

A Newton–Krylov method with an approximate analytical Jacobian for implicit solution of Navier–Stokes equations on staggered overset-curvilinear grids with immersed boundaries

The explicit and semi-implicit schemes in flow simulations involving complex geometries and moving boundaries suffer from time-step size restriction and low convergence rates. Implicit schemes can be used to overcome these restrictions, but implementing them to solve the Navier–Stokes equations is not straightforward due to their non-linearity. Among the implicit schemes for non-linear equations, Newton-based techniques are preferred over fixed-point techniques because of their high convergence rate but each Newton iteration is more expensive than a fixed-point iteration. Krylov subspace methods are one of the most advanced iterative methods that can be combined with Newton methods, i.e., Newton–Krylov Methods (NKMs) to solve non-linear systems of equations. The success of NKMs vastly depends on the scheme for forming the Jacobian, e.g., automatic differentiation is very expensive, and matrix-free methods without a preconditioner slow down as the mesh is refined. A novel, computationally inexpensive analytical Jacobian for NKM is developed to solve unsteady incompressible Navier–Stokes momentum equations on staggered overset-curvilinear grids with immersed boundaries. Moreover, the analytical Jacobian is used to form a preconditioner for matrix-free method in order to improve its performance. The NKM with the analytical Jacobian was validated and verified against Taylor–Green vortex, inline oscillations of a cylinder in a fluid initially at rest, and pulsatile flow in a 90 degree bend. The capability of the method in handling complex geometries with multiple overset grids and immersed boundaries is shown by simulating an intracranial aneurysm. It was shown that the NKM with an analytical Jacobian is 1.17 to 14.77 times faster than the fixed-point Runge–Kutta method, and 1.74 to 152.3 times (excluding an intensively stretched grid) faster than automatic differentiation depending on the grid (size) and the flow problem. In addition, it was shown that using only the diagonal of the Jacobian further improves the performance by 42–74% compared to the full Jacobian. The NKM with an analytical Jacobian showed better performance than the fixed point Runge–Kutta because it converged with higher time steps and in approximately 30% less iterations even when the grid was stretched and the Reynold number was increased. In fact, stretching the grid decreased the performance of all methods, but the fixed-point Runge–Kutta performance decreased 4.57 and 2.26 times more than NKM with a diagonal and full Jacobian, respectivley, when the stretching factor was increased. The NKM with a diagonal analytical Jacobian and matrix-free method with an analytical preconditioner are the fastest methods and the superiority of one to another depends on the flow problem. Furthermore, the implemented methods are fully parallelized with parallel efficiency of 80–90% on the problems tested. The NKM with the analytical Jacobian can guide building preconditioners for other techniques to improve their performance in the future.

Analysis of a Decoupled Time-Stepping Scheme for Evolutionary Micropolar Fluid Flows

Micropolar fluid model consists of Navier-Stokes equations and microrotational velocity equations describing the dynamics of flows in which microstructure of fluid is important. In this paper, we propose and analyze a decoupled time-stepping algorithm for the evolutionary micropolar flow. The proposed method requires solving only one uncoupled Navier-Stokes and one microrotation subphysics problem per time step. We derive optimal order error estimates in suitable norms without assuming any stability condition or time step size restriction.

A Resistive Boundary Condition Enhanced DGTD Scheme for the Transient Analysis of Graphene

In this paper, the electromagnetic (EM) features of graphene are characterized by a discontinuous Galerkin time-domain (DGTD) algorithm with a resistive boundary condition (RBC). The atomically thick graphene is equivalently modeled using an RBC by regarding the graphene as an infinitesimally thin conductive sheet. To incorporate RBC into the DGTD analysis, the surface conductivity of the graphene composed of contributions from both intraband and interband terms is first approximated by rational basis functions using the fast-relaxation vector-fitting (FRVF) method in the Laplace domain. Next, through the inverse Laplace transform, the corresponding time-domain matrix equations in integral can be obtained. Finally, these matrix equations are solved by time-domain finite integral technique (FIT). For elements not touching the graphene sheet, however, the well-known Runge-Kutta (RK) method is employed to solve the two first-order time-derivative Maxwell’s equations. The application of the surface boundary condition significantly alleviates the memory consuming and the limitation of time step size required by Courant–Friedrichs–Lewy (CFL) condition. To validate the proposed algorithm, various numerical examples are presented and compared with available references.

Fast Transient Analysis of Nonuniform Multiconductor Transmission Lines Using HIE-Block-LIM

This letter describes a hybrid implicit-explicit block latency insertion method (HIE-block-LIM) for the fast simulation of nonuniform multiconductor transmission lines (MTLs). In the HIE-block-LIM, an implicit difference method is used with respect to the current variables in the y-direction, and an explicit method is adopted to update the other variables. The HIE-block-LIM can alleviate a time step size limitation of the existing block-LIM by taking both advantages of the explicit and implicit difference methods. Numerical results show that the HIE-block-LIM is suitable for the fast simulation of the nonuniform MTLs.

Locally Implicit Time Integration Strategies in a Discontinuous Galerkin Method for Maxwell’s Equations

An attractive feature of discontinuous Galerkin (DG) spatial discretization is the possibility of using locally refined space grids to handle geometrical details. However, locally refined meshes lead to severe stability constraints on explicit integration methods to numerically solve a time-dependent partial differential equation. If the region of refinement is small relative to the computational domain, the time step size restriction can be overcome by blending an implicit and an explicit scheme where only the solution variables living at fine elements are treated implicitly. The downside of this approach is having to solve a linear system per time step. But due to the assumed small region of refinement relative to the computational domain, the overhead will also be small while the solution can be advanced in time with step sizes determined by the coarse elements. In this paper, we present two locally implicit time integration methods for solving the time-domain Maxwell equations spatially discretized with a DG method. Numerical experiments for two-dimensional problems illustrate the theory and the usefulness of the implicit–explicit approaches in presence of local refinements.

A Kinetic Model for Gas Mixtures Based on a Fokker-Planck Equation

We present a generalized nonlinear Fokker–Planck equation, which describes the dynamics of rarefied monatomic gas mixture flows. The devised kinetic model leads to correct transfer of energy and momentum between gas species and to consistent evolution of molecular stresses and heat fluxes with respect to the generalized Boltzmann equation. Thus, the correct diffusion coefficient together with the mixture viscosity and mixture heat conductivity coefficients are obtained. The strength of the presented model lies on the computational efficiency, which is due to the fact that the resulting stochastic processes are continuous in time. Therefore, unlike in Direct Simulation Monte Carlo (DSMC), here simulated particles do not collide with each other, but move along independent continuous stochastic paths. Another aspect of the new Fokker-Planck model equation is that here the effect of collisions is described via drift and diffusion type processes. Accordingly, a scheme can be derived for which the time step size limitation of the corresponding numerical simulation becomes independent of the Knudsen number. Consequently, this leads to more efficient simulations, especially in low or intermediate Knudsen numbers. Results are presented for helium-argon mixture in a one dimensional geometry. The calculated mixture viscosity is found to be in accordance with experimental data, which reveals the accuracy and relevance of the approach.

Preprocess static subdomain decomposition in practical cases of 2D unsteady hydraulic simulation

Explicit finite volume methods are frequently used and widely accepted in hydraulic models based on the shallow water approximation. The main drawback of the approach is the time step size limit imposed by the Courant–Friedrichs–Lewy numerical stability constraint. This leads to excessively long computational times in large scale cases of practical interest. At the same time, the accuracy of the numerical results is associated to the use of fine computational meshes able to achieve enough spatial resolution. Taking into account that hydraulic modelers do not have access, in general, to large computational facilities, suitable and useful parallelization techniques are required. Furthermore, if high performance computing facilities are used, it is usually necessary to provide an estimation of the requirements of computational load to cover the length of the simulation. In this work the suitability of a preprocess static subdomain decomposition is explored and presented as a promising strategy to improve the efficiency of 2D unsteady shallow water computational models over dry bed in medium scale computational facilities and, at the same time, is useful to provide a preprocess computational time estimation if large scale computational facilities are going to be used.

Dynamic Subtiming-Based Implicit Nonoscillating Scheme for Contaminant Transport Modeling

A dynamic subtime-stepping method is described for solving contaminant transport problems that utilize higher-order implicit time-marching procedures with higher-accuracy nonoscillating spatial-discretization methods to resolve sharp-plume fronts in advection-dominated systems. Nonoscillating spatial-discretization methods for the advective term prevent unphysical oscillations and minimize numerical diffusion. Second-order temporal accuracy is achieved by using the Crank-Nicholson implicit scheme, however nonoscillating properties may be violated for Courant Numbers larger than one, leading to spurious oscillations. The proposed subtiming method allows use of small time-step sizes in critical portions of a domain, with larger time-step sizes in other locations. This locally limits the Courant number where required and still keeps the general solution free of time-step size restrictions. This technique makes it possible to apply higher-order nonoscillating schemes with higher-order temporal weighting for advection-dominated flows, even when the Courant number is much greater than one. Feasibility and applicability of the dynamic implicit subtiming method are demonstrated through three proof-of-concept example problems.

The discrete equation method (DEM) for fully compressible, two-phase flows in ducts of spatially varying cross-section

For the simulation of light water nuclear reactor coolant flows, general two-phase models (valid for all volume fractions) have been generally used which, while allowing for velocity disequilibrium, normally force pressure equilibrium between the phases (see, for example, the numerous models of this type described in (Städtke, 2006)). These equations are not hyperbolic, their physical wave dynamics are incorrect, and their solution algorithms rely on dubious truncation error induced artificial viscosity to render them numerically well-posed over a portion of the computational spectrum. The inherent problems of the traditional approach to multiphase modeling, which begins with an averaged system of (ill-posed) partial differential equations (PDEs) which are then discretized to form a numerical scheme, are avoided by employing a new homogenization method known as the discrete equation method (DEM) (Abgrall and Saurel, 2003). This method results in well-posed hyperbolic systems, this property being important for transient flows. This also allows a clear treatment of non-conservative terms (terms involving interfacial variables and volume fraction gradients) permitting the solution of interface problems without conservation errors, this feature being important for the direct numerical simulation of two-phase flows.Unlike conventional methods, the averaged system of PDEs for the mixture are not used, and the DEM method directly obtains a well-posed discrete equation system from the single-phase conservation laws, producing a numerical scheme which accurately computes fluxes for arbitrary number of phases. The method effectively uses a sequence of single-phase Riemann problem solutions. Phase interactions are accounted for by Riemann solvers at each interface. Non-conservative derivative product terms are correctly approximated. Some of the closure relations missing from the traditional approach are automatically obtained. Lastly, the continuous equation system resulting from the discrete equations can be identified by taking the continuous limit with weak-wave assumptions. In this work, this approach is tested by constructing a DEM model for the flow of two compressible phases in 1-D ducts of spatially varying cross-section with explicit time integration. An analytical equation of state is included for both water vapor and liquid phases, and a realistic interphase mass transfer model is developed based on interphase heat transfer. A robust compliment of boundary conditions is developed and discussed. Though originally conceived as a first step toward implict time integration of the DEM method (to relieve time step size restrictions due to stiffness and to achieve tighter coupling of equations) in multidimensions, this model offers some unique capabilities for incorporation into next generation light water reactor safety analysis codes. We demonstrate, on a converging-diverging two-phase nozzle, that this well-posed, 2-pressure, 2-velocity DEM model can be integrated to a realistic and meaningful steady state with both phases treated as compressible.