Predictor-corrector Integration Research Articles

Multi-step Hermite-Birkhoff predictor-corrector schemes

Multi-step Hermite-Birkhoff predictor-corrector schemes

Modelling deformable boundary by spherical particle for normal contact

The normal contact of the elastic spherical particle with deformable boundary is investigated in terms of the Discrete Element Method (DEM). The particle of the prescribed radius is moving under gravity and the initial velocity. The deformable boundary is treated as rigidly fixed spherical particle with variable elasticity modulus and variable radius. The limit case, approaching the infinite radius presents an elastic half-space, while increasing of the elasticity modulus presents the rigid boundary, respectively.
 The linear model and the nonlinear Hertz contact model used in the discrete element method are investigated numerically by applying the 5th-order Gear’s predictorcorrector integration scheme. The numerical model is tested by comparing it with analytical solution. The time variations of the particle positions, velocities and accelerations are presented. On the basis of simulation results the limit values of the boundary particle parameters are evaluated and recommendations for the boundary article parameters required in DEM simulation are drown.

Inertial instability and phase error in Euler forward predictor-corrector time integration schemes: Improvement of modeling Great Lakes thermal structure and circulation using FVCOM

This study investigates the inertial stability properties and phase error of numerical time integration schemes in several widely-used ocean and atmospheric models. These schemes include the most widely used centered differencing (i.e., leapfrog scheme or the 3-time step scheme at n-1, n, n+1) and 2-time step (n, n+1) 1st-order Euler forward schemes, as well as 2nd-stage and 3rd- and 4th-stage Euler predictor-corrector (PC) schemes. Previous work has proved that the leapfrog scheme is neutrally stable with respect to the Coriolis force, with perfect inertial motion preservation, an amplification factor (AF) equal to unity, and a minor overestimation of the phase speed. The 1st-order Euler forward scheme, on the other hand, is known to be unconditionally inertially unstable since its AF is always greater than unity. In this study, it is shown that 3rd- and 4th-order predictor-corrector schemes 1) are inertially stable with weak damping if the Coriolis terms are equally split to n+1 (new value) and n (old value); and 2) introduce an artificial computational mode. The inevitable phase error associated with the Coriolis parameter is analyzed in depth for all numerical schemes. Some schemes (leapfrog and 2nd-stage PC schemes) overestimate the phase speed, while the others (1st-order Euler forward, 3rd- and 4th-stage PC schemes) underestimate it. To preserve phase speed as best as possible in a numerical model, alternating a scheme that overestimates the phase speed with a scheme that underestimates the phase speed is recommended. Considering all properties investigated, the leapfrog scheme is still highly recommended for a time integration scheme. As an example, a comparison between a leapfrog scheme and a 1st-order Euler forward scheme is presented to show that the leapfrog scheme reproduces much better vertical thermal stratification and circulation in the weakly-stratified Great Lakes.

Step forward on nonlinear differential equations with the Atangana–Baleanu derivative: Inequalities, existence, uniqueness and method

Step forward on nonlinear differential equations with the Atangana–Baleanu derivative: Inequalities, existence, uniqueness and method

Discrete exterior calculus discretization of two-phase incompressible Navier-Stokes equations with a conservative phase field method

We present a discrete exterior calculus (DEC) based discretization scheme for incompressible two-phase flows. Our physically-compatible exterior calculus discretization of single phase flow is extended to simulate immiscible two-phase flows with discontinuous changes in fluid properties such as density and viscosity across the interface. The two-phase incompressible Navier-Stokes equations and conservative phase field equation for interface capturing are first transformed into the exterior calculus framework. The discrete counter part of these smooth equations is obtained by substituting with discrete differential forms and discrete exterior operators. We prove the boundedness of the method for the first order Euler forward and predictor-corrector time integration schemes in the DEC framework. With a proper choice of two free parameters, the scheme remains phase field bounded without the requirement of any ad hoc mass redistribution. We verify the scheme against several standard test cases (for interface capturing) comprising not only the flat domains but also the curved domains, leveraging the advantage that DEC operators are independent of the coordinate system. The results show excellent properties of boundedness, mass conservation and convergence. Moreover, we demonstrate the ability of the scheme towards handling large density and viscosity ratios as well as surface tension in the simulation of various two phase flow physical phenomena on flat or curved surfaces.

SAFARI – A low and hyperthermal energy ion scattering simulation

SAFARI [1][2] is a classical molecular dynamics simulation designed to assist with analysis of low and hyperthermal energy ion spectra (1 eV - 1 keV) for ion-surface scattering experiments. SAFARI is optimized for this specific energy range to reduce the required runtime, and to produce more precise spectra. This paper discusses the algorithms used in SAFARI, as well as various optimizations implemented for the specific ranges involved. Program Title: SAFARI CPC Library link to program files: https://doi.org/10.17632/t5kxgnmjwv.1 Developer's repository link: https://github.com/SINS-Lab/SEA-SAFARI Code Ocean capsule: https://codeocean.com/capsule/3787536 Licensing provisions: MIT Programming language: C++ Nature of problem: SAFARI models low and hyperthermal energy (1 eV - 1 keV) ion scattering spectra. This energy range requires molecular dynamics techniques to adequately reproduce experimental spectra, as the energies are too low for appropriate use of the sequential binary collision approximation. Several approximations and optimisations are used to reduce the runtime of the simulation, to allow generating larger statistical datasets. Solution method: SAFARI utilises a predictor-corrector integrator, and classical Hamiltonian mechanics to model the trajectory of the ions. Ion-lattice forces are computed by a combination of image parameters, optional frictional losses, and binary pair potentials. The pair potentials are only computed for nearby lattice sites to reduce the computational complexity of the simulation.

Discrete Exterior Calculus Discretization of Two-Phase Incompressible Navier-Stokes Equations with a Conservative Phase Field Method

We present a discrete exterior calculus (DEC) based discretization scheme for incompressible two-phase flows. Our physically-compatible exterior calculus discretization of single phase flow is extended to simulate immiscible two-phase flows with discontinuous changes in fluid properties such as density and viscosity across the interface. The two-phase incompressible Navier-Stokes equations and conservative phase field equation for interface capturing are first transformed into the exterior calculus framework. The discrete counter part of these smooth equations is obtained by substituting with discrete differential forms and discrete exterior operators. We prove the boundedness of the method for the first order Euler forward and predictor-corrector time integration schemes in the DEC framework. With a proper choice of two free parameters, the scheme remains phase field bounded without the requirement of any ad hoc mass redistribution. We verify the scheme against several standard test cases (for interface capturing) comprising not only the flat domains but also the curved domains, leveraging the advantage that DEC operators are independent of the coordinate system. The results show excellent properties of boundedness, mass conservation and convergence. Moreover, we demonstrate the ability of the scheme towards handling large density and viscosity ratios as well as surface tension in the simulation of various two phase flow physical phenomena on flat or curved surfaces.

Ab Initio Direct Dynamics.

The reactivity and dynamics of molecular systems can be explored computationally by classical trajectory calculations. The traditional approach involves fitting a functional form of a potential energy surface (PES) to the energies from a large number of electronic structure calculations and then integrating numerous trajectories on this fitted PES to model the molecular dynamics. The ever-decreasing cost of computing and continuing advances in computational chemistry software have made it possible to use electronic structure calculations directly in molecular dynamics simulations without first having to construct a fitted PES. In this "on-the-fly" approach, every time the energy and its derivatives are needed for the integration of the equations of motion, they are obtained directly from quantum chemical calculations. This approach started to become practical in the mid-1990s as a result of increased availability of inexpensive computer resources and improved computational chemistry software. The application of direct dynamics calculations has grown rapidly over the last 25 years and would require a lengthy review article. The present Account is limited to some of our contributions to methods development and various applications. To improve the efficiency of direct dynamics calculations, we developed a Hessian-based predictor-corrector algorithm for integrating classical trajectories. Hessian updating made this even more efficient. This approach was also used to improve algorithms for following the steepest descent reaction paths. For larger molecular systems, we developed an extended Lagrangian approach in which the electronic structure is propagated along with the molecular structure. Strong field chemistry is a rapidly growing area, and to improve the accuracy of molecular dynamics in intense laser fields, we included the time-varying electric field in a novel predictor-corrector trajectory integration algorithm. Since intense laser fields can excite and ionize molecules, we extended our studies to include electron dynamics. Specifically, we developed code for time-dependent configuration interaction electron dynamics to simulate strong field ionization by intense laser pulses. Our initial application of ab initio direct dynamics in 1994 was to CH2O → H2 + CO; the calculated vibrational distributions in the products were in very good agreement with experiment. In the intervening years, we have used direct dynamics to explore energy partitioning in various dissociation reactions, unimolecular dissociations yielding three fragments, reactions with branching after the transition state, nonstatistical dynamics of chemically activated molecules, dynamics of molecular fragmentation by intense infrared laser pulses, selective activation of specific dissociation channels by aligned intense infrared laser fields, angular dependence of strong field ionization, and simulation of sequential double ionization.

Explicit meshfree {{varvec{u}}}-{{varvec{p}}}_mathbf{mathrm{w}} solution of the dynamic Biot formulation at large strain

In this paper, an efficient and robust methodology to simulate saturated soils subjected to low-medium frequency dynamic loadings under large deformation regime is presented. The coupling between solid and fluid phases is solved through the dynamic reduced formulation u-p_mathrm{w} (solid displacement – pore water pressure) of the Biot’s equations. The additional novelty lies in the employment of an explicit two-steps Newmark predictor-corrector time integration scheme that enables accurate solutions of related geomechanical problems at large strain without the usually high computational cost associated with the implicit counterparts. Shape functions based on the elegant Local Maximum Entropy approach, through the Optimal Transportation Meshfree framework, are considered to solve numerically different dynamic problems in fluid saturated porous media.

Macroscopic constitutive model for ergodic and non-ergodic lead-free relaxors

A fully electromechanically coupled, three dimensional phenomenological constitutive model for relaxor ferroelectric materials was developed for the use in a finite-element-method (FEM) solution procedure. This macroscopic model was used to simulate the macroscopic electromechanical response of lead-free ergodic [Formula: see text] and non-ergodic [Formula: see text] relaxor materials. The presented constitutive model is capable of accounting for the observed pinched hysteretic response as well as non-deviatoric polarization induced strain and internal order transitions. Time integration of the history dependent internal variables is done with a predictor-corrector integration scheme. The adaptability of the constitutive model regarding the pinching of the hystereses is shown. Simulations are compared to experimental observations.

Re-integration with anchor points algorithm for ab initio molecular dynamics.

A new integration scheme for ab initio molecular dynamics (MD) is proposed in this work for efficient propagation using large time steps (e.g., 2.0fs or a larger time step with one ab initio evaluation of gradients for the dynamics point and one additional evaluation for the anchor point per dynamics step). This algorithm is called re-integration with anchor points (REAP) integrator. The REAP integrator starts from a quadratic potential energy surface based on the updated Hessian to propagate the system to the halfway of the MD step that is called the anchor point. Then, an approximate dynamics position for this step is obtained by the propagation based on an interpolated surface using the anchor point and the previous MD point. The approximate dynamics step can be further improved by the re-integration steps, i.e., integration based on the interpolated surface using the calculated energies, gradients, and updated Hessians of the previous step, the anchor point, and the approximate current step. A trajectory only needs one analytical Hessian calculation at the initial geometry, and thereafter, only calculations of gradients are required. This integrator can be considered either as a generalization of Hessian-based predictor-corrector integration with substantial improvement of accuracy and efficiency or as a dynamics on interpolated surfaces that are built on the fly. An automatic correction scheme is implemented by comparing the interpolated energies and gradients to the actual ones to ensure the quality of the interpolations at a certain level. The tests in this work show that the REAP method can increase computational efficiency by more than one order of magnitude than that of the velocity Verlet integrator and more than twice that of Hessian-based predictor-corrector integration.

Radial integral boundary element method for simulating phase change problem with mushy zone

A radial integral boundary element method (BEM) is used to simulate the phase change problem with a mushy zone in this paper. Three phases, including the solid phase, the liquid phase, and the mushy zone, are considered in the phase change problem. First, according to the continuity conditions of temperature and its gradient on the liquid-mushy interface, the mushy zone and the liquid phase in the simulation can be considered as a whole part, namely, the non-solid phase, and the change of latent heat is approximated by heat source which is dependent on temperature. Then, the precise integration BEM is used to obtain the differential equations in the solid phase zone and the non-solid phase zone, respectively. Moreover, an iterative predictor-corrector precise integration method (PIM) is needed to solve the differential equations and obtain the temperature field and the heat flux on the boundary. According to an energy balance equation and the velocity of the interface between the solid phase and the mushy zone, the front-tracking method is used to track the move of the interface. The interface between the liquid phase and the mushy zone is obtained by interpolation of the temperature field. Finally, four numerical examples are provided to assess the performance of the proposed numerical method.

Nonlinear free vibration analysis of piezoelectric laminated plate with random actuation electric potential difference and thermal loading

• Solutions to moderately large amplitude vibrations of piezoelectric thin laminated plates. • Mechanical-thermal-electrical coupled system of nonlinear equations. • Displacement-stress mixed formulation. • Random actuation electric potential and thermal load on nonlinear frequency ratio. Solutions to nonlinear free vibration analysis of an arbitrarily laminated thin rectangular plate subjected to thermo-piezoelectric load is investigated. The governing equation of motion is based on von Kármán’s deflection of thin plates to capture the effects of moderate geometric nonlinearity . Induced strain actuation assumption is confined to a linear model for applied low electric fields and varying temperature distribution along the thickness direction. The solution for lateral displacement field for simply supported and clamped condition edges are determined by displacement-stress; mixed formulation. The applied Galerkin’s projection reduces the governing nonlinear partial differential equation to a time dependent ordinary differential equation . Solutions are obtained using a predictor-corrector numerical integration scheme. The effects of randomly applied actuation potential and thermal loading on nonlinear frequency ratio are examined. The marginal values provided using various probabilistic methods are found to be in good agreement with each other.

Creep damage and fracture of notched specimens under static and fast periodic loading

Creep-damage processes in notched specimens subjected to static and periodic loading are studied experimentally. Subsequent simulations were carried out using Rabotnov evolution equation. In case of dynamic creep (fast periodic loading) the constitutive equations derived by use of asymptotic and averaging methods were used. Numerical results were obtained by the combination of FEM and predictor-corrector time integration scheme. Fracture process was studied using designed numerical procedure based on the elimination of the finite elements with critical values of damage parameter. The fracture times and directions of the crack propagation were determined. The qualitative difference between processes of crack propagation in cases of static and periodic loading were observed.

Real-time time-dependent density functional theory using density fitting and the continuous fast multipole method.

An implementation of real-time time-dependent density functional theory (RT-TDDFT) within the TURBOMOLE program package is reported using Gaussian-type orbitals as basis functions, second and fourth order Magnus propagator, and the self-consistent field as well as the predictor-corrector time integration schemes. The Coulomb contribution to the Kohn-Sham matrix is calculated combining density fitting approximation and the continuous fast multipole method. Performance of the implementation is benchmarked for molecular systems with different sizes and dimensionalities. For linear alkane chains, the wall time for density matrix time propagation step is comparable to the Kohn-Sham (KS) matrix construction. However, for larger two- and three-dimensional molecules, with up to about 5,000 basis functions, the computational effort of RT-TDDFT calculations is dominated by the KS matrix evaluation. In addition, the maximum time step is evaluated using a set of small molecules of different polarities. The photoabsorption spectra of several molecular systems calculated using RT-TDDFT are compared to those obtained using linear response time-dependent density functional theory and coupled cluster methods.

How fast do young star clusters expel their natal gas? Estimating the upper limit of the gas expulsion time-scale

ABSTRACT Formation of massive stars within embedded star clusters starts a complex interplay between their feedback, inflowing gas, and stellar dynamics, which often includes close stellar encounters. Hydrodynamical simulations usually resort to substantial simplifications to model embedded clusters. Here, we address the simplification which approximates the whole star cluster by a single sink particle, which completely neglects the internal stellar dynamics. In order to model the internal stellar dynamics, we implement a Hermite predictor–corrector integration scheme to the hydrodynamic code flash. As we illustrate by a suite of tests, this integrator significantly outperforms the current leap-frog scheme, and it is able to follow the dynamics of small compact stellar systems without the necessity to soften the gravitational potential. We find that resolving individual massive stars instead of representing the whole cluster by a single energetic source has a profound influence on the gas component: for clusters of mass less than $\approx3 \times 10^3 \, \mathrm{M}_{\odot }$ , it slows gas expulsion by a factor of ≈5 to $\approx 1 \, \mathrm{Myr}$, and it results in substantially more complex gas structures. With increasing cluster mass (up to $\approx 3\times 10^3 \, \mathrm{M}_{\odot }$), the gas expulsion time-scale slightly decreases. However, more massive clusters ($\gtrsim 5\times 10^3 \, \mathrm{M}_{\odot }$) are unable to clear their natal gas with photoionizing radiation and stellar winds only if they form with a star formation efficiency (SFE) of 1/3. This implies that the more massive clusters are either cleared with another feedback mechanism or they form with an SFE higher than 1/3.

Nonlinear Dynamics of Swinging Clapper Bells under Arbitrary or Resonant Forcing Functions

Study of swinging clapper bells involves aspects encompassing sound and acoustic engineering, mechanical engineering, and structural engineering. From the musical point of view, clapper bells are directly played idiophone instruments, where the playing device, the clapper, although directly excited, is not explicitly controlled by the bell ringer. The achievement of a clear and optimal sound mainly depends on the acoustic characteristics of the bell and on the regularity of the clapper strokes, which is not only governed by the ringing style and the relevant parameters of clapper and bell but also by the real time corrections to the excitation introduced by trained bell ringers. In fact, despite centuries of experience allowed to optimize the bell performances, standardizing proportions and mounting arrangements, effective sound control requires some fine tuning of the forcing function. Another crucial topic, especially in view of assessing existing structures, regards the evaluation of time histories of the actions transmitted by the bell to the pivots and the study of the interactions between the bell and the supporting structures, belfries, and bell-towers. “Ringability” of swinging bells and bell-structure interactions are usually tackled in the framework of rigid body dynamics, so arriving at an initial value problem, governed by a system of two second order nonlinear ordinary differential equations (ODEs), whose solutions are piecewise-defined functions. In the relevant literature, numerical solutions of the system are commonly sought using built-in algorithms provided in advanced software packages; since the use of such general algorithms is subject to some restrictions, especially regarding the forcing functions, validity of the results is often limited. The present study focuses on an innovative procedure to solve the equations of motion. The method, extremely fast and effective, is based on original numerical explicit-implicit predictor-corrector integration algorithms with constant time step, duly validated reproducing the outcomes of relevant reference case studies. Each time the clapper strikes the bell a new “piece” of the solution is initialized, so avoiding user interventions in the elaboration phase. Independently on the oscillation amplitude and on the duration of the considered time interval, the algorithms can successfully manage undamped oscillations; friction and viscosity damped oscillations; free oscillations in transient and stationary phases; and can be applied also to solve stiff equations. Furthermore, the capability of the proposed methods to deal with arbitrary forcing functions is particularly innovative. The outcomes of relevant case studies, regarding the oscillations of the old tenor bell of the Great St. Mary church in Cambridge, confirm the potentialities of the method, also highlighting some topical issues, involving, for example, the assessment of damping equivalence. Finally, a pioneering feature of the algorithms is their ability to handle and to define “resonant” forcing functions, continuously tuning the frequency of the excitation to the natural frequency of the oscillation, according to the oscillation amplitude.

Squeeze film dampers supporting high-speed rotors: Rotordynamics

This work develops a finite element based multi-mass flexible rotor model for theoretical investigation of the influence of the squeeze film damper lubricant inertia on the unbalance-induced steady-state and transient vibration amplitudes of high speed turbomachinery. The rotordynamic model is developed by applying the principles of finite element analysis to discretize the rotor components, including the rotor shaft and disk, into local elements with mass, stiffness, and gyroscopic matrices. Subsequently, the local matrices are assembled together to develop the global model of the rotordynamic system. The influence of squeeze film damper lubricant inertia is incorporated into the model by using short-length cavitated damper models with retaining springs executing circular-centered orbits. Additionally, the rotordynamic model incorporating the nonlinear squeeze film damper models is iteratively solved in the time domain by applying a predictor-corrector transient modal integration numerical method and the steady-state and transient motions of the rotor system are investigated under different rotor and squeeze film damper parameters. The results of the study verify the substantial influence of squeeze film damper lubricant inertia on attenuating the vibrations of high-speed turbomachinery. Furthermore, the developed rotordynamic model delivers an efficient and powerful platform for the analysis of high-speed turbomachinery, including jet engines and gas turbines.

Three-step Predictor-Corrector Finite Element Schemes for Consolidation Equation

An accurate, stable, and efficient three-step predictor-corrector time integration method is considered, for the first time, to obtain numerical solution for the one-dimensional consolidation equation within a finite and spectral element framework. Theoretical order of accuracy and stability conditions are provided. The three-step predictor-corrector time integration method is third-order accurate and shows a larger stability region than the forward Euler method when applied to the one-dimensional consolidation equation. Furthermore, numerical results are in agreement with analytical solutions previously derived by the authors.

A Comprehensive Framework for Coupled Nonlinear Aeroelasticity and Flight Dynamics of Highly Flexible Aircrafts

A framework to model and analyze the coupled nonlinear aeroelasticity and flight dynamics of highly flexible aircrafts is presented. The methodology is based on the dynamics of 3D co-rotational beams. The coupling of axial, bending and torsional effects is added to the stiffness and mass matrices of Euler–Bernoulli beam to capture the most relevant characteristics of a real wing structure. The finite-state aerodynamic model is coupled with the structural model to simulate the unsteady aerodynamics. A scheme of mixed end-point and mid-point time-marching algorithms is proposed and applied into the implicit predictor–corrector integration, where the end-point algorithm is used in the predictor step for efficiency and mid-point algorithm in corrector step for accuracy. The ground, body and airflow axes for flight dynamics are re-defined by the global and elemental ones for structural dynamics, followed by the redefinitions of local Euler angles and airflow angles of each element. The framework can be used for quick analyses of flexible aircrafts in conceptual and preliminary design phases, including linear and nonlinear trim, aerodynamic load estimation, stability assessment, time-domain simulations and flight performance evaluations. The results show the payload mass and its distributions will significantly affect the trim state and longitudinal stability of highly flexible aircrafts.