Time Integration Method Research Articles

Explicit Runge–Kutta Numerical Manifold Method for Solving the Burgers’ Equation via the Hopf–Cole Transformation

This paper presents an efficient numerical manifold method for solving the Burgers’ equation. To improve accuracy and streamline the solution process, we apply a nonlinear function transformation technique that reformulates the original problem into a linear diffusion equation. We utilize a dual cover mesh along with an explicit multi-step time integration method for spatial and temporal discretization, respectively. Constant cover functions are employed across mathematical covers, interconnected by a linear weight function for each manifold element. The full discretization formulation is derived using the Galerkin weak form. To efficiently compute the inverse of the symmetric positive definite mass matrix, we employ the Crout algorithm. The performance and convergence of our method are rigorously evaluated through several benchmark numerical tests. Extensive comparisons with exact solutions and alternative methods demonstrate that our approach delivers an accurate, stable, and efficient computational scheme for the Burgers’ equation.

A fully probabilistic framework to compute the residual rockfall risk in presence of mitigation measures

AbstractRockfall events are expected to rise throughout the future due to climate change and extreme meteorological events. In the perspective of climate change adaptation, an accurate quantification of the risk is needed, together with a precise assessment of the effectiveness of protective measures eventually installed. All the possible block detachment scenarios together with their occurrence probability should be considered, and a time span should be selected. A fully probabilistic framework to compute the risk in absence and in presence of a protective structure is herein proposed, and a time-integrated reliability-based method, developed by the authors, is applied to define the failure probability of the protective measure. The complete method, in absence and presence of a rockfall barrier, is applied to a study case, and the residual risk in presence of the barrier is quantified. The results show the importance of considering all the possible detachment situations to have reliable results in terms of both risk and effectiveness of the protective measure quantification.

Acoustic and soliton propagation using fully-discrete energy preserving partially implicit scheme in homogeneous and heterogeneous mediums

This study presents an energy preserving partially implicit scheme for the simulation of wave propagation in homogeneous and heterogeneous mediums. Despite its implicit nature, the developed scheme does not require any explicit numerical or analytical inversion of the coefficient matrix. Theoretical analysis and numerical experiments are performed to validate the energy preserving properties of the fully-discrete scheme. Convergence analysis is also performed to assess the rate of convergence of the developed scheme. The efficiency and accuracy of the developed scheme are validated by numerical solutions of wave propagation in layered heterogeneous mediums. Furthermore, simulations of soliton propagation following nonlinear sine-Gordon and Klein-Gordon equations in homogeneous and heterogeneous mediums are discussed. Numerical solutions are also compared with the results available in the literature. The present method accurately resolves the physical characteristics of the chosen problems, competing well with existing multi-stage time-integration methods. Moreover, it has significantly lower computational complexity than the four-stage, fourth-order Runge-Kutta-Nyström method.

An elasto-visco-plastic constitutive model for snow: Theory and finite element implementation

Snow exhibits unique mechanical behaviour due to its evolving properties influenced by temperature, stress conditions, and viscous effects. This paper introduces a nonlinear constitutive model for snow, featuring new formulations for the yield function and strain rate potential, and incorporating viscosity, sintering, and degradation effects. The model is numerically integrated into the Abaqus/Standard FEM software using a fully implicit backward Euler method for time integration and Powell’s hybrid method for linearizing the nonlinear system of equations. The robustness and stability of the numerical scheme ensure accurate simulation of snow behaviour under various loading conditions. The model is finally validated against experimental data available in the literature, demonstrating its effectiveness and reliability in capturing the complex mechanical response of snow.

An explicit time integration method for Boussinesq approximation

AbstractThis paper presents a novel approach to the discretization of the Boussinesq equation using finite elements, coupled with an explicit pressure correction scheme. The entire numerical scheme is formulated exclusively in terms of matrix‐vector products, enabling a highly efficient numerical solution. The key advantage of this approach lies in its suitability for fast computations and seamless parallelization on graphics processing units (GPUs). We present the numerical scheme and discuss different computational examples.

Numerical study of bifurcation structures in reflected shock-wave/laminar-boundary-layer interaction within an end-wall tube

Abstract The primary aim of this study is to analyze the unsteady characteristics of the interaction between a reflected shock wave and a laminar boundary layer in an end-wall shock tube. Our direct numerical simulations at shock Mach numbers of Ms=1.9, 2.5, and 3.5 using a fifth-order WENO scheme and three-step Runge-Kutta time integration method revealed inhomogeneity and anisotropy in the shock bifurcation. Surprisingly, the upper and lower bifurcated structures maintain a notably asymmetric flow during the forward propagation of the reflected shock bifurcation. The inverse flow in the bifurcation resembles a crooked earthworm structure, exhibiting high-frequency oscillations indicative of instability. However, at higher shock intensities, the earthworm transforms into a stable strip-like configuration, facilitating the entrapment of inverse flow and leading to rapid bifurcation height growth and early convergence. Additionally, isolated islands with high density, temperature, and pressure emerge in the transitional region behind the bifurcated shocks, due to variations in wave propagation speed.

Spatially mixed implicit–explicit schemes in hydro-mechanically coupled soil dynamics

Mixed implicit–explicit (IMEX) schemes combine elements of both implicit and explicit time integration methods. This allows to efficiently handle boundary value problems with severe differences in stiffness, such as often encountered in soil dynamics. A spatially mixed IMEX scheme, in which different parts of the computational domain may be treated implicitly or explicitly depending on their behaviour, is discussed in this work. Domain decomposition by means of a mortar contact discretisation is used for this purpose and multiple systems of equations are solved. The merits of the scheme in solving problems of soil dynamics are demonstrated, for which the governing equations are extended to hydro-mechanical coupling.

A novel iterative algorithm for nonlinear inelastic dynamic analysis of truss structures; MCB-Spline+Sp time integration method

The focus of this work is on the development of a novel algorithm for the nonlinear dynamic analysis of structures. In the sake of a general and easy-to-implement solution method, the tangential stiffness matrix is calculated by combining the Cauchy-Riemann equations derived from the differentiation of complex functions and using the nodal displacement perturbation vector. An implicit direct time integration method based on a cubic B-Spline is integrated with the incremental-iterative method using a quadrature rule based on the interpolation of spline functions, which is referred to as the MCB-Spline+Sp time integration method. A step-by-step algorithm for developing a computer code to implement this method is provided. The effectiveness of the proposed approach is evaluated through several nonlinear time history analyses, which demonstrate its accuracy and efficiency by observing less computational efforts and fast convergence rate. Moreover, a comparison is made between the developed method and well-established techniques such as the Newmark and Standard-Bathe time integration methods in combination with the iterative Newton-Raphson method.

DG-IMEX method for a two-moment model for radiation transport in the [formula omitted] limit

We consider neutral particle systems described by moments of a phase-space density and propose a realizability-preserving numerical method to evolve a spectral two-moment model for particles interacting with a background fluid moving with nonrelativistic velocities. The system of nonlinear moment equations, with special relativistic corrections to O(v/c), expresses a balance between phase-space advection and collisions and includes velocity-dependent terms that account for spatial advection, Doppler shift, and angular aberration. The model is conservative for the correct O(v/c) Eulerian-frame number density and is consistent, to O(v/c), with Eulerian-frame energy and momentum conservation. This model is closely related to the one promoted by Lowrie et al. [1] and similar to models currently used to study transport phenomena in large-scale simulations of astrophysical environments. The proposed numerical method is designed to preserve moment realizability, which guarantees that the moments correspond to a nonnegative phase-space density. The realizability-preserving scheme consists of the following key components: (i) a strong stability-preserving implicit-explicit (IMEX) time-integration method; (ii) a discontinuous Galerkin (DG) phase-space discretization with carefully constructed numerical fluxes; (iii) a realizability-preserving implicit collision update; and (iv) a realizability-enforcing limiter. In time integration, nonlinearity of the moment model necessitates solution of nonlinear equations, which we formulate as fixed-point problems and solve with tailored iterative solvers that preserve moment realizability with guaranteed global convergence. We also analyze the simultaneous Eulerian-frame number and energy conservation properties of the semi-discrete DG scheme and propose a “spectral redistribution” scheme that promotes Eulerian-frame energy conservation. Through numerical experiments, we demonstrate the accuracy and robustness of this DG-IMEX method and investigate its Eulerian-frame energy conservation properties.

Implicit Peer Triplets in Gradient-Based Solution Algorithms for ODE Constrained Optimal Control

It is common practice to apply gradient-based optimization algorithms to numerically solve large-scale ODE constrained optimal control problems. Gradients of the objective function are most efficiently computed by approximate adjoint variables. High accuracy with moderate computing time can be achieved by such time integration methods that satisfy a sufficiently large number of adjoint order conditions and supply gradients with higher orders of consistency. In this paper, we upgrade our former implicit two-step Peer triplets constructed in [Algorithms, 15:310, 2022] to meet those new requirements. Since Peer methods use several stages of the same high stage order, a decisive advantage is their lack of order reduction as for semi-discretized PDE problems with boundary control. Additional order conditions for the control and certain positivity requirements now intensify the demands on the Peer triplet. We discuss the construction of 4-stage methods with order pairs (3, 3) and (4, 3) in detail and provide three Peer triplets of practical interest. We prove convergence of order s-1\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$s-1$$\\end{document}, at least, for s-stage methods if state, adjoint and control satisfy the corresponding order conditions. Numerical tests show the expected order of convergence for the new Peer triplets.

Hydrodynamics of multicomponent vesicles: A phase-field approach

In this paper, we introduce a thermodynamically-consistent phase-field model to investigate the hydrodynamics of inextensible multicomponent vesicles in various fluid flows with inertial forces. Our model couples the fluid field, surface concentration field representing chemical species on the membrane, and vesicle dynamics, while enforcing global area and volume constraints through a Lagrange multiplier method. Specifically, we employ full Navier–Stokes equations for the fluid field, the Cahn–Hilliard equations for the species concentration on the membrane, a nonlinear advection–diffusion equation to describe the evolution of the vesicle membrane, and an additional equation to enforce local inextensibility. We utilize a residual-based variational multiscale method for the Navier–Stokes equations and a standard Galerkin finite element framework for the remaining equations. The PDEs are solved using an implicit, monolithic scheme based on the generalized-α time integration method. We extend previous models for homogeneous vesicles (Aland et al., 2014; Valizadeh and Rabczuk, 2022) to multicomponent vesicles, introducing Cahn–Hilliard equations while maintaining thermodynamic consistency. Additionally, we employ isogeometric analysis (IGA) for higher accuracy. We present a variety of two-dimensional numerical examples, including multicomponent vesicles in shear flow and Poiseuille flow with and without obstructions, using the resistive immersed surface method to handle obstructions. Furthermore, we provide three-dimensional simulations of multicomponent vesicles in Poiseuille flow.

Non-smooth self-excited vibration of a novel dynamical model for a disc brake

This paper proposes a new dynamic model for the study of friction-induced self-excited vibration of a disc brake system, where the pad’s motions in both radial and circumferential/tangential directions are included, which is in stark contrast to the previous studies that normally consider the pad’s motion in the tangential/circumferential direction only. The non-smooth dynamics of the system including three different states of motion, i.e., stick, slip and separation, is investigated. Both the linear stability analysis and the transient dynamic analysis are performed. The numerical results in the linear stability analysis indicate that the inclusion of pad’s radial motion in the present brake model significantly expands the ranges of operating parameters for dynamic instability than the brake model with only circumferential/tangential motion for the pad. For the transient dynamic analysis, two different methods, i.e., the time integration method and the shooting method, are employed for the calculation of steady-state response. The accuracy and efficiency of the shooting method are subsequently examined. The numerical results show rich bifurcation behaviours of the steady-state response in the present brake model with the variations of brake pressure and disc speed , and (the stiffness of the inclined spring in the radial direction) is a key parameter for controlling the occurrence of chaotic vibration in the system. Conflict of Interest Statement The authors declare no conflicts of interest.

On the measurement of the nonlinear dynamics of sandwich sector plate surrounded by the auxetic concrete foundation: Introducing a machine learning algorithm for nonlinear problems

The measurement of the nonlinear dynamics of sandwich rotary sector plates is crucial for measurement engineers as it enables the precise analysis and understanding of the behavior of advanced composite structures under dynamic conditions. These measurements help in identifying and mitigating potential issues related to vibration, stability, and fatigue, which are critical in industries like aerospace, automotive, and civil engineering. Accurate data on nonlinear dynamics aids in the optimization of design, enhancing performance, durability, and safety. Furthermore, it supports the development of predictive maintenance strategies, reducing downtime and operational costs, and contributing to the advancement of engineering materials and techniques. So, in the current work, for the first time, nonlinear dynamics of sandwich rotary sector plate via a mathematical modeling and machine learning algorithm is presented. First, due to a lack of information on the nonlinear dynamics of sandwich rotary sector plates, a dataset for training, testing, and validating the machine learning algorithm is presented. For this purpose, using Hamilton’s principle, Von-Karman nonlinearity, and first-order shear deformation theory (FSDT), the nonlinear governing equations are obtained. Also, due to increasing stiffness and stability, an auxetic concrete foundation covers the sandwich structure. After that using the two-dimensional generalized differential quadrature method (GDQM) and Newmark’s time integration method (NTIM), the nonlinear equations are solved. This research provides valuable insights for engineers in designing advanced composite structures with improved dynamic properties. The results also contribute to the broader understanding of nonlinear dynamic interactions in complex material systems, paving the way for innovative applications in various engineering fields.

Forced vibration characteristics of viscoelastic variable stiffness laminated composite plates using time and frequency domain approaches

The linear forced vibration characteristics of viscoelastic variable stiffness laminated composite plates (VVSLC) are studied using a generalized Maxwell model and finite element method. The integral form of the viscoelastic constitutive relation is converted to the incremental form for finite element formulation based on the Reissner–Mindlin plate theory. The recursive relations are developed to compute the current time-step solution using only the previous time-step solution. The periodic response directly in the time domain is obtained using shooting technique coupled with Newmark’s time integration method. The implementation of shooting technique for Boltzmann integral-based viscoelasticity for curvilinear fibre composite plates is done for the first time in this study. For the comparison purpose, the response/resonance frequency/modal loss factor is also obtained using equivalent complex modulus based viscoelastic correspondence principle. It is observed that the variation in fibre orientation and boundary conditions leads to significant variations in response, stress/moment resultant amplitude and the damping factor of the VVSLC plates. Further, the present time domain based approach is capable of predicting damping factor at all forcing frequencies whereas complex eigenvalue analysis can predict damping factor only at discrete resonance frequency. Based on the detailed studies, it is found that the curvilinear fibre composite plates depict a significant reduction in response/moment resultants compared to straight fibre composite plates.

A dual-modified implicit time integration method for three-dimensional impact modelling within the framework of the SBFEM

This paper presents a dual-modified implicit time integration method (DMITIM) for three-dimensional impact modelling within the framework of the scaled boundary finite element method (SBFEM). Controlled by a single parameter, this method effectively dissipates high-frequency responses in impact problems without affecting low-frequency responses, ensuring time integration stability. To address the issue of false oscillations during impacts, additional Lagrange multipliers are introduced to supplement the velocity and acceleration corrections at each time integration step. For nonmatching meshes, by leveraging the mesh flexibility of the SBFEM, a rapid and straightforward remeshing method is utilized to establish a node-to-node contact scheme. The numerical results demonstrate that the DMITIM significantly reduces false oscillations in the velocity and contact force during impact events. It also showcases the flexibility of the SBFEM in handling nonmatching meshes. The combination of the SBFEM with the DMITIM provides enhanced simulation accuracy and stability for impact problems compared to the results from ABAQUS using the HHT-α time integration algorithm.

On Properties of Adjoint Systems for Evolutionary PDEs

We investigate the geometric structure of adjoint systems associated with evolutionary partial differential equations at the fully continuous, semi-discrete, and fully discrete levels and the relations between these levels. We show that the adjoint system associated with an evolutionary partial differential equation has an infinite-dimensional Hamiltonian structure, which is useful for connecting the fully continuous, semi-discrete, and fully discrete levels. We subsequently address the question of discretize-then-optimize versus optimize-then-discrete for both semi-discretization and time integration, by characterizing the commutativity of discretize-then-optimize methods versus optimize-then-discretize methods uniquely in terms of an adjoint-variational quadratic conservation law. For Galerkin semi-discretizations and one-step time integration methods in particular, we explicitly construct these commuting methods by using structure-preserving discretization techniques.

(Battery Student Slam 8 Award Winner) A Comparison of Different Convergence Criteria to Improve the Performance of a P2D Modelling Algorithm of a Lithium-Ion Battery

Lithium-ion batteries (LiBs) are used as the storage for electric vehicles, grid energy and consumer electronics mainly due to the LiBs’ high energy density, long cycle life, and reliable manufacturing. Effective performance and diagnostic management of LiBs is important and is implemented via Battery Management Systems (BMS) [1]. A BMS should also maximise the safety and lifetime of battery packs and therefore, needs to estimate the internal battery states including the state of charge (SoC), state of health (SoH), state of power (SoP), state of function (SoF), and state of safety (SoS) [2]. Accurate estimation of the battery SoX can be realized with the help of a battery model; an algorithm representing the electrochemical, mechanical, electrical, and thermal behaviour of the LiB cells.Cell-level battery models typically fall into three categories: (i) Empirical Models, (ii) Equivalent Circuit Models (ECM) and (iii) Physics-based Electrochemical Models (EM). Electrochemical Models represent the internal cell processes within a cell based on its governing fundamental physical and electrochemical equations [3].The Partial two-dimensional (P2D) model is considered a high-fidelity Electrochemical Model (EM), which in its most basic form, captures three major physio-chemical processes in a cell: mass transfer via reduction/oxidisation (redox) reactions, diffusion of neutral lithium in the electrodes and diffusion of lithium-ions in the electrolyte. The P2D model captures the different physio-chemical phenomena using a set of coupled Partial Differential Algebraic Equations (PDAEs) [4].Despite accuracy, P2D models are complex mathematical models with typically no analytical solution [5]. Therefore, to implement the full-order version of the P2D model, numerical-iterative solver algorithms are required. Multiple works have looked at different methods of optimizing and solving the P2D model. These methods include using different discretisation schemes, such as the finite difference method and the finite volume method, using different iterative solvers and different convergence criteria [5]. However, the different works on iterative solvers we have reviewed to date [5], haven’t investigated the effects of the different criteria on improving the performance of P2D iterative solvers.Hence, in this work, we aim to build on our previous work [6], where we proposed an algorithm to implement the P2D model using the finite volume spatial discretisation method (FVM) and Euler time integration method. We then analyse how the role of (i) the choice of convergence criteria, (ii) the choice of initial guess and (iii) the choice of root-finding method used to converge to the convergence criteria, can be varied to improve the performance of the iterative solver.Therefore, the contributions of this work can be summarised as: An algorithm that implements the isothermal P2D model of a LiB using the Finite Volume Method and Euler integration methods.An analysis of the changes in algorithm performance when the following performance factors are varied; (i) The convergence criteria, (ii) the choice of initial guess and (iii) the choice of root-finding method.An optimal configuration of performance criteria that solves the isothermal-P2D model for a 1C discharge simulation, with an average speed of 3 seconds and a root mean square error of less than 1% when compared to the commercial solver COMSOL. [1] C. P. Grey and D. S. Hall, “Prospects for lithium-ion batteries and beyond—a 2030 vision,” Nature Communications 2020 11:1, vol. 11, no. 1, pp. 1–4, Dec. 2020, doi: 10.1038/s41467-020-19991-4.[2] L. Lu, X. Han, J. Li, J. Hua, and M. Ouyang, “A review on the key issues for lithium-ion battery management in electric vehicles,” Journal of Power Sources, vol. 226. pp. 272–288, Mar. 15, 2013. doi: 10.1016/j.jpowsour.2012.10.060.[3] S. Abada, G. Marlair, A. Lecocq, M. Petit, V. Sauvant-Moynot, and F. Huet, “Safety focused modeling of lithium-ion batteries: A review,” J Power Sources, vol. 306, pp. 178–192, Feb. 2016, doi: 10.1016/J.JPOWSOUR.2015.11.100.[4] M. Doyle, T. F. Fuller, and J. Newman, “Modeling of Galvanostatic Charge and Discharge of the Lithium/Polymer/Insertion Cell,” J Electrochem Soc, vol. 140, no. 6, pp. 1526–1533, Jun. 1993, doi: 10.1149/1.2221597/XML.[5] A. M. Ramos, “On the well-posedness of a mathematical model for Lithium-ion batteries,” Appl Math Model, vol. 40, no. 1, pp. 115–125, May 2015, doi: 10.48550/arxiv.1506.00605.[6] T. Wickramanayake, M. Javadipour, and K. Mehran, “A Novel Root-Finding Algorithm to Solve the Pseudo-2D Model of a Lithium-ion Battery,” 2023 IEEE International Conference on Electrical Systems for Aircraft, Railway, Ship Propulsion and Road Vehicles and International Transportation Electrification Conference, ESARS-ITEC 2023, 2023, doi: 10.1109/ESARS-ITEC57127.2023.10114840.

A Semi-Explicit Algorithm for Parameters Estimation in a Time-Fractional Dual-Phase-Lag Heat Conduction Model

This paper presents a new semi-explicit algorithm for parameters estimation in a time-fractional generalization of a dual-phase-lag heat conduction model with Caputo fractional derivatives. It is shown that this model can be derived from a general linear constitutive relation for the heat transfer by conduction when the heat conduction relaxation kernel contains the Mittag–Leffler function. The model can be used to describe heat conduction phenomena in a material with power-law memory. The proposed algorithm of parameters estimation is based on the time integral characteristics method. The explicit representations of the thermal diffusivity and the fractional analogues of the thermal relaxation time and the thermal retardation are obtained via a Laplace transform of the temperature field and utilized in the algorithm. An implicit relation is derived for the order of fractional differentiation. In the algorithm, this relation is resolved numerically. An example illustrates the proposed technique.

New implicit time integration schemes combining high frequency damping with high second order accuracy

It is shown that weighted linear combinations of the generalised-α method and certain related higher order schemes allow for the formulation of unconditionally stable single step time integration methods with improved second order accuracy and more targeted high-frequency damping. It is also shown that, if the user controlled high frequency damping parameter ρ∞ is set to zero, the new schemes can be expressed as linear multistep backward difference formulae and, in a particular case, recover Park's method. The performance of the proposed methods is illustrated in terms of mathematical analysis and a number of linear and nonlinear numerical examples including finite element based solutions of the incompressible Navier-Stokes equations.

High order numerical treatment of the Riesz space fractional advection-dispersion equation via matrix transform method

In this paper, a mixed high order finite difference scheme-Padé approximation method is applied to obtain numerical solution of the Riesz space fractional advection-dispersion equation(RSFADE). This method is based on the high order finite difference scheme that derived from fractional centered difference and Padé approximation method for space and time integration, respectively. The stability analysis of the proposed method is discussed via theoretical matrix analysis. Numerical experiments are presented to confirm the theoretical results of the proposed method.