Explicit Numerical Integration Research Articles

Towards modelling AR Sco: generalized particle dynamics and strong radiation-reaction regimes

ABSTRACT Numerical simulations of relativistic plasmas have become more feasible, popular, and crucial for various astrophysical sources with the availability of computational resources. The necessity for high-accuracy particle dynamics is especially highlighted in pulsar modelling due to the extreme associated electromagnetic fields and particle Lorentz factors. Including the radiation-reaction force in the particle dynamics adds even more complexity to the problem, but is crucial for such extreme astrophysical sources. We have also realized the need for such modelling concerning magnetic mirroring and particle injection models proposed for AR Sco, the first white dwarf pulsar. This paper demonstrates the benefits of using higher-order explicit numerical integrators with adaptive time-step methods to solve the full particle dynamics with radiation-reaction forces included. We show that for standard test scenarios, namely various combinations of uniform E- and B-fields and a static dipole B-field, the schemes we use are equivalent to and in extreme field cases outperform standard symplectic integrators in accuracy. We show that the higher-order schemes have massive computational time improvements due to the adaptive time-steps we implement, especially in non-uniform field scenarios and included radiation reaction where the particle gyro-radius rapidly changes. When balancing accuracy and computational time, we identified the adaptive Dormand–Prince eighth-order scheme to be ideal for our use cases. The schemes we use maintain accuracy and stability in describing the particle dynamics and we indicate how a charged particle enters radiation-reaction equilibrium and conforms to the analytical Aristotelian Electrodynamics expectations.

Phase-field description of fracture in NiTi single crystals

A phase-field model for thermomechanically-induced fracture in NiTi at the single crystal level, i.e., fracture under loading paths that may take advantage of either of the functional properties of NiTi–superelasticity or shape memory effect–, is presented, formulated within the kinematically linear regime. The model accounts for reversible phase transformation from austenite to martensite habit plane variants and plastic deformation in the austenite phase. Transformation-induced plastic deformation is viewed as a mechanism for accommodation of the local deformation incompatibility at the austenite–martensite interfaces and is accounted for by introducing an interaction term in the free energy derived based on the Mori–Tanaka and Kröner micromechanical assumptions and the hypothesis of martensite instantaneous growth within austenite. Based on experimental observations suggesting that NiTi fractures in a stress-controlled manner, damage is assumed to be driven by the elastic energy, i.e., phase transformation and plastic deformation are assumed to contribute in crack formation and growth indirectly through stress redistribution. The model is restricted to quasistatic mechanical loading (no latent heat effects), thermal loading sufficiently slow with respect to the time rate of heat transfer by conduction (no thermal gradients), and a temperature range below Md, which is the temperature above which the austenite phase is stable, i.e., stress-induced martensitic transformation is suppressed. The numerical implementation of the model is based on an efficient scheme of viscous regularization in both phase transformation and plastic deformation, an explicit numerical integration via a tangent modulus method, and a staggered scheme for the coupling of the unknown fields. The model is shown able to capture transformation-induced toughening, i.e., stable crack advance attributed to the shielding effect of inelastic deformation left in the wake of the growing crack under nominal isothermal loading, actuation-induced fracture under a constant bias load, and crystallographic dependence on crack pattern.

Numerical analysis of non-proportional biaxial reverse experiments with a two-surface anisotropic cyclic plasticity-damage approach

This paper deals with the numerical analysis of ductile damage and fracture behavior under non-proportional biaxial reverse loading conditions. A two-surface anisotropic cyclic elastic–plastic-damage continuum model is adequately presented, which takes into account the Bauschinger effect, the stress-differential effect, and the change of hardening rate after reverse loading. An efficient Euler explicit numerical integration algorithm, based on the inelastic (plastic or plastic-damage) predictor-elastic corrector approach, is utilized to analyze the stress and finite strain loading histories. Detailed discussions are provided on different numerical integration-related consistent tangent operators that achieve convergence within the global Newton–Raphson scheme. The proposed continuum model is implemented into the commercial software Ansys as a user-defined subroutine (UMAT). Furthermore, the novel non-proportional biaxial tensile reverse experiments are performed to validate the proposed continuum model. The associated numerical simulations investigate the stability and accuracy of the proposed algorithm and material model.

On stress and deformation accumulation of continuously welded rails under cyclic thermal loading in high-speed railways

This paper focuses on the stress and deformation accumulation of continuously welded rails (CWR) in ballastless tracks under cyclic thermal loading, which could pose a potential threat to the stability and service life of CWR but related theoretical studies are yet to be conducted. Firstly, the nonlinear longitudinal resistance of fasteners (LRF) resulting from the CWR-rail pad contact interface is mathematically described by extending the classical Dhal friction model and verified via an in-door experiment. Both the hysteretic and softening properties of the LRF can be captured compared with the traditional ideal elastoplastic model. Most importantly, we originally propose a quasi-dynamic method (QDM) to tackle nonlinear static problems by extending an explicit numerical integration algorithm from the dynamics field to the statics. Afterwards, the quasi-dynamic equilibrium equations of CWR in the tunnel transition zone and on the simply-supported beam bridges subjected to thermal loads and LRF are constructed, and the reliability of QDM to obtain static responses is well-validated by existing analytical methods. Finally, the limitation of the ideal elastoplastic LRF model, distribution characteristics and accumulation process of rail longitudinal deformation and axial additional stress are fully unveiled. The current work is conducive to better recognition of the longitudinal mechanical behavior of CWR under cyclic thermal loads.

Splitting Methods for Semi-Classical Hamiltonian Dynamics of Charge Transfer in Nonlinear Lattices

We propose two classes of symplecticity-preserving symmetric splitting methods for semi-classical Hamiltonian dynamics of charge transfer by intrinsic localized modes in nonlinear crystal lattice models. We consider, without loss of generality, one-dimensional crystal lattice models described by classical Hamiltonian dynamics, whereas the charge (electron or hole) is modeled as a quantum particle within the tight-binding approximation. Canonical Hamiltonian equations for coupled lattice-charge dynamics are derived, and a linear analysis of linearized equations with the derivation of the dispersion relations is performed. Structure-preserving splitting methods are constructed by splitting the total Hamiltonian into the sum of Hamiltonians, for which the individual dynamics can be solved exactly. Symmetric methods are obtained with the Strang splitting of exact, symplectic flow maps leading to explicit second-order numerical integrators. Splitting methods that are symplectic and conserve exactly the charge probability are also proposed. Conveniently, they require only one solution of a linear system of equations per time step. The developed methods are computationally efficient and preserve the structure; therefore, they provide new means for qualitative numerical analysis and long-time simulations for charge transfer by nonlinear lattice excitations. The properties of the developed methods are explored and demonstrated numerically considering charge transport by mobile discrete breathers in an example model previously proposed for a layered crystal.

New Analytical Results and Comparison of 14 Numerical Schemes for the Diffusion Equation with Space-Dependent Diffusion Coefficient

We examine the one-dimensional transient diffusion equation with a space-dependent diffusion coefficient. Such equations can be derived from the Fokker–Planck equation and are essential for understanding the diffusion mechanisms, e.g., in carbon nanotubes. First, we construct new, nontrivial analytical solutions with the classical self-similar Ansatz in one space dimension. Then we apply 14 different explicit numerical time integration methods, most of which are recently introduced unconditionally stable schemes, to reproduce the analytical solution. The test results show that the best algorithms, especially the leapfrog-hopscotch, are very efficient and severely outperform the conventional Runge–Kutta methods. Our results may attract attention in the community who develops multi-physics engineering software.

FREQUENCY FEATURES OF THE NUMERICAL METHOD SAMPLING OF DIGITAL CONTROL SYSTEMS

Context. The studies of the frequency properties of the explicit multistep numerical integrators which use for sampling of continuous transfer function in the digital control systems, are conducted in this article. Numerical integrators in such systems implement as an integral parts of the digital regulators. Objective. The goal of this research is an analysis of the behavior of explicit numerical integrators of different orders, which areused to discretize continuous systems, in order to study their impact on the properties of the synthesized digital system. Method. Numerical methods of integration are considered as digital filters, the behavior of which is studied by the frequency characteristics method. To do this, the z-transform apparatus was used. Integrators’ discrete transfer functions were found for frequency analysis using the Control Systems Toolbox package of the mathematical application MATLAB. For further analysis, two closed feedback test structures were used: with integrators in the forward channel and in the feedback loop. Both variants of structures were studied by the frequency characteristics method for sampling using numerical integrators of 1st–6th orders. Results. The inefficiency of using high-order numerical integrators for continuous systems’ discretization is shown. Given the behavior of the frequency characteristics of test systems, the most rational is the use of low-order integrators, namely – the first and second orders. Establishing the cause of this phenomenon requires additional research, in particular, to identify the possible impact of additional zeros and poles of discrete transfer functions of the numerical integrators. Conclusions. The use of low-order integrators, namely the first and second orders, is the most rational for sampling of digital control systems and the inefficiency of using high-order numerical integrators to sample continuous systems is proven.

Modelling unsaturated soils with the Material Point Method. A discussion of the state-of-the-art

Many natural hazards involve large deformations of unsaturated soils, e.g. rainfall-induced landslides, embankment collapses due to wetting, seepage-induced instabilities of dams and levees, etc. The study of these phenomena requires accounting for the complex hydro-mechanical interactions between solid skeleton and pore fluids and modelling large deformations to predict the post-failure behaviour, which poses significant computational challenges. In recent years, several hydro-mechanical coupled MPM formulations were developed to model saturated and unsaturated soils. These approaches are slightly different in terms of governing equations, integration schemes and have been implemented in different MPM software; thus, they benefit from various computational features. The purpose of this paper is to present an overview of the available MPM approaches to model unsaturated soils discussing differences and similarities of the formulations and their impact on the results under different conditions in a range of geotechnical applications. In addition, the effect of partially saturated conditions on the critical time step in explicit numerical integration schemes is studied for the first time. Different analytical expressions are derived and compared with the numerical results.

On the Applicability of Approximate Rolling and Sliding Contact Algorithms in Anisothermal Problems with Thermal Softening

The residual stress approximation methods formulated by McDowell and Moyar, Jiang and Sehitoglu, and McDowell for rolling and sliding contact problems are reconsidered in the context of single anisothermal loading cycles and isotropic hardening. A consistent extention to incorporate thermal softening is developed and the generalized thermoelastoplastic algorithms are cast into a proper predictor–corrector formulation. Detailed explicit and implicit numerical integration strategies are presented and validated using specifically designed finite element models that conform to the underlying mechanical assumptions. Then, the applicability of the approximate algorithms to anisothermal problems with isotropic hardening and thermal softening is analyzed by assuming a rate-independent Johnson–Cook-type yield stress model and by comparing the obtained transient and residual stresses to results from full-scale finite element half-space models under varying loading and strain-hardening intensities. An in-depth, comparative discussion on the adequacy of the algorithms in conjunction with the justification of their respective mechanical simplifications follows. Sufficiently strong strain hardening is found to be a prerequisite for accurate predictions, and Jiang and Sehitoglu’s approach is deemed to be preferable for the considered type of problem. The conclusions drawn from the investigations are discussed in the context of common applications with particular emphasis on manufacturing process modeling and the corresponding guidelines are proposed for such use cases.

Comparison of implicit and explicit numerical integration schemes for a bounding surface soil model without elastic range

Standard integration schemes for rate constitutive equations were designed within the classical theory of plasticity. Consequently, they rely on the assumption that a yield criterion defines a range of purely elastic material behaviour. Many constitutive models for non-cohesive soils discard that assumption. They account for inelastic deformations without employing a yield criterion. This simplifies the formulation of the models but raises questions concerning their numerical implementation. In particular, it is not fully clear which algorithmic method is most appropriate for the numerical integration of constitutive models without elastic range. To investigate this, two different stress point algorithms for a critical state bounding surface model for sands were developed and implemented. The explicit method employs substepping to automatically control the local error. The implicit method updates stresses and state variables through a local Newton iteration, the Jacobian of which is computed by numerical differentiation. The two algorithms were compared by means of calculations at integration point level as well as with respect to a boundary value problem. The results show that for a given level of accuracy, the explicit update procedure is significantly more efficient than the implicit one. This holds regardless of initial state parameters and input strain increment magnitudes.

A Numerical Study for Solving Mathematical Fuzzy Population Model

Population models are used to determine maximum harvest for agriculturists and to understand the dynamics of biological invasions and for environmental conservation. Also, population models are used to understand the spread of parasites, viruses and disease. However, the mathematical model of population can be used to make future prediction of the growth of population. In this paper, numerical study for solving mathematical fuzzy population model has been introduced. The development of direct explicit numerical integrators of RK-type to be consistent with solving fuzzy ordinary differential equations. The computational efficiency of the numerical method has been developed. The advantage and the efficiency of the proposed method have been shown clearly using the numerical results which are agree well with analytical solutions due to the proposed method is more efficient and accurate method.

Developed RKM Method for Solving Ninth-Order Ordinary Differential Equations with Applications

In this study, the development of direct explicit numerical integrators of RK-type for solving a class of ninth-order ordinary differential equations (ODEs) to develop the computational efficiency of the methods. The objective of this paper is to generalize the numerical integrators of RK type for solving classes of ODEs of orders up to eighth-order ODEs. Using the approach of Taylor-expansion, we have derived the order-conditions (OCs) for RKM integrators. Based on these OCs, developed method of ninth-order with five-stage has been derived. The zero stability of RKM integrator is proven. Stability polynomial of RKM integrator for test problem has been introduced. The advantage and the efficiency of the proposed method have been shown clearly using the numerical results which are agree well with analytical solutions due to the proposed method is more efficient and accurate method.

Structure-Preserving Numerical Integrators for Hodgkin--Huxley-Type Systems

Motivated by the Hodgkin--Huxley model of neuronal dynamics, we study explicit numerical integrators for “conditionally linear” systems of ordinary differential equations. We show that splitting and composition methods, when applied to the Van der Pol oscillator and to the Hodgkin--Huxley model, do a better job of preserving limit cycles of these systems for large time steps, compared with the “Euler-type” methods (including Euler's method, exponential Euler, and semi-implicit Euler) commonly used in computational neuroscience, with no increase in computational cost. These limit cycles are important to preserve, due to their role in neuronal spiking. Splitting methods even compare favorably to the explicit exponential midpoint method, which is twice as expensive per step. The second-order Strang splitting method is seen to perform especially well across a range of nonstiff and stiff dynamics.

GeCo: Geometric Conservative nonstandard schemes for biochemical systems

We generalize the nonstandard Euler and Heun schemes in order to provide explicit geometric numerical integrators for biochemical systems, here denoted as GeCo schemes, that preserve both positivity of the solutions and linear invariants. We relax the request on the order convergence of the denominator function for the first-order approximation and we let it depend on the step size also throughout the solution approximating values. The first-order variant is exact on a two-dimensional linear test problem. Moreover, we introduce a class of modified mGeCo(α) schemes and, by tuning the parameter α≥1, we improve the numerical performance of GeCo integrators on some examples taken from the literature. A numerical comparison with BBKS and mBBKS schemes, which are implicit integrators, positive and conserve linear invariants, show the gain in efficiency of both GeCo and mGeCo(α) procedures as they generate similar errors with an explicit functional form.

Numerical integrators based on the Magnus expansion for nonlinear dynamical systems

Explicit numerical integration algorithms up to order four based on the Magnus expansion for nonlinear non-autonomous ordinary differential equations are presented and tested on problems possessing qualitative (very often, geometric) features that is convenient to preserve under numerical discretization. The range of applications covers augmented dynamical systems, highly oscillatory problems and nonlinear non-autonomous partial differential equations of evolution previously discretized in space.

Explicit, fully implicit and forward gradient numerical integration of a hyperelasto-viscoplastic constitutive model for amorphous polymers undergoing finite deformation

Following the growing use of amorphous polymers in an expanding range of applications, interest for numerical modeling of polymer behavior has greatly increased. Together with reliable constitutive models, stable, accurate and rapid integration algorithms valid for large deformations need to be developed. Here, in the framework of hyperelasto-viscoplasticity and multiplicative split formulation, three integration algorithms (explicit, fully implicit and forward gradient) are generated for a constitutive polymer model and respective stability is investigated. The algorithms are furthermore implemented in a commercial Finite Element code and simulation of a full field tensile test is shown to capture the actual deformation behavior of polymers.

Symmetric, explicit numerical integrator for molecular dynamics equations of motion with a generalized friction

A general mathematical scheme to construct symmetric, explicit numerical integrators of Newtonian equations of motion endowed with a generalized friction is provided for a molecular dynamics (MD) study. The exact integrations are done for all the decomposed vector fields, including the one that contains the friction term. On the basis of the symmetric composition scheme with the adjoint for the resulting maps, integrators with any local order of accuracy can be systematically constructed. Among them, the second order P2S1 integrator gives the least evaluation of atomic force and potential, which are most time consuming in MD simulations. As examples of the friction function, three functional types are considered: constant, Laurent polynomial, and exponential with respect to the kinetic energy. Several MD equations of motion fall into these categories, and the numerical examinations of their integrators using model systems give positive results on the accuracy and efficiency. The extended phase-space scheme also presents an invariant function, which allows us to easily detect numerical errors in the integration process by monitoring the function value.

Critical time step for DEM simulations of dynamic systems using a Hertzian contact model

SummaryThe discrete element method (DEM) typically uses an explicit numerical integration scheme to solve the equations of motion. However, like all explicit schemes, the scheme is only conditionally stable, with the stability determined by the size of the time step. Currently, there are no comprehensive techniques for estimating appropriate DEM time steps when a nonlinear contact interaction is used. It is common practice to apply a large factor of safety to these estimates to ensure stability, which unnecessarily increases the computational cost of these simulations. This work introduces an alternative framework for selecting a stable time step for nonlinear contact laws, specifically for the Hertz‐Mindlin contact law. This approach uses the fact that the discretised equations of motion take the form of a nonlinear map and can be analysed as such. Using this framework, we analyse the effects of both system damping and the initial relative velocity of collision on the critical time step for a Hertz‐Mindlin contact event between spherical particles.

Analysis of Numerical Methods to Include Dynamic Constraints in an Optimal Power Flow Model

The optimization of the operation of power systems including steady state and dynamic constraints is efficiently solved by Transient Stability Constrained Optimal Power Flow (TSCOPF) models. TSCOPF studies extend well-known optimal power flow models by introducing the electromechanical oscillations of synchronous machines. One of the main approaches in TSCOPF studies includes the discretized differential equations that represent the dynamics of the system in the optimization model. This paper analyzes the impact of different implicit and explicit numerical integration methods on the solution of a TSCOPF model and the effect of the integration time step. In particular, it studies the effect on the power dispatch, the total cost of generation, the accuracy of the calculation of electromechanical oscillations between machines, the size of the optimization problem and the computational time.

Simulation of spatial multibody systems with friction

ABSTRACTA formulation for modeling and simulation of friction effects in spatial multibody systems is presented. Constraint reaction forces on rigid bodies that are connected by joints that support friction are derived as functions of Lagrange multipliers, using D’Alembert’s principle. Friction forces acting on bodies are calculated as a function of joint geometry, constraint reaction forces that are functions of Lagrange multipliers, and relative velocities at constraint contact points that are determined by system kinematics. Friction forces are implemented in index 0 differential-algebraic equations of motion that are solved numerically using explicit and implicit numerical integration methods. Spatial examples are presented, yielding accurate results and demonstrating that the systems are not stiff, even in the presence of friction and stiction.