Numerical Time Integration Research Articles

A meshless fragile points method for the solution of the monodomain model for cardiac electrophysiology simulation

Meshless methods for in silico modeling and simulation of cardiac electrophysiology are gaining more and more popularity. These methods do not require a mesh and are more suitable than the Finite Element Method (FEM) to simulate the activity of complex geometrical structures like the human heart. However, challenges such as numerical integration accuracy and time efficiency remain, which limit their applicability. Recently, the Fragile Points Method (FPM) has been introduced in the meshless methods family. It uses local, simple, polynomial, discontinuous functions to construct trial and test functions in the Galerkin weak form. This allows for accurate integration and improved efficiency while enabling the imposition of essential and natural boundary conditions as in FEM. In this work, we consider the application of FPM for cardiac electrophysiology simulation. We derive the cardiac monodomain model using the FPM formulation and we solve several benchmark problems in 2D and 3D. We show that FPM leads to solutions of accuracy and efficiency similar to FEM while alleviating the need for a mesh. Additionally, FPM demonstrates comparable convergence to FEM in the considered benchmarks.

Bounds for rotating convection at infinite Prandtl number from semidefinite programs

The technique of semidefinite programming is used to replace direct numerical simulations and time integrations of turbulent flows with an optimization problem which is easier to solve but which returns only upper bounds to the time averages the direct numerical simulation would have produced. This idea is applied to rotating Rayleigh-B\'enard convection at infinite Prandtl number. The obtained bounds are not sharp but they reproduce the qualitative behavior of the results known from direct simulations and experiments.

A Multirate Discontinuous-Galerkin-in-Time Framework for Interface-Coupled Problems

A framework is presented to design multirate time stepping algorithms for two dissipative models with coupling across a physical interface. The coupling takes the form of boundary conditions imposed on the interface, relating the solution variables for both models to each other. The multirate aspect arises when numerical time integration is performed with different time step sizes for the component models. In this paper, we seek to identify a unified approach to develop multirate algorithms for these coupled problems. This effort is pursued though the use of discontinuous-Galerkin time stepping methods, acting as a general unified framework, with different time step sizes. The subproblems are coupled across user-defined intervals of time, called coupling windows, using polynomials that are continuous on the window. The coupling method is shown to reproduce the correct interfacial energy dissipation, discrete conservation of fluxes, and asymptotic accuracy. In principle, methods of arbitrary order are possible. As a first step, herein we focus on the presentation and analysis of monolithic methods for advection-diffusion models coupled via generalized Robin-type conditions. The monolithic methods could be computed using a Schur-complement approach. The framework and analysis herein accommodate multirate coupling for a wide variety of single-step, multistep, and multistage methods. A concrete example of Crank--Nicolson integrators coupled in a multirate fashion is also included to illustrate the ideas. We conclude with some discussion of future developments, such as different interface conditions and partitioned methods.

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.

GPU-accelerated time integration of Gross-Pitaevskii equation with discrete exterior calculus

The quantized vortices in superfluids are modeled by the Gross-Pitaevskii equation whose numerical time integration is instrumental in the physics studies of such systems. In this paper, we present a reliable numerical method and its efficient GPU-accelerated implementation for the time integration of the three-dimensional Gross-Pitaevskii equation. The method is based on discrete exterior calculus which allows us the usage of more versatile spatial discretization than traditional finite difference and spectral methods are applicable to. We discretize the problem using six different natural crystal structures and observe the correct choices of spatial tiling to decrease the truncation error and increase the reliability compared to Cartesian grids. We pay attention to the computational performance optimizations of the GPU implementation and measure speedups of up to 152-fold when compared to a reference CPU implementation. We parallelize the implementation further to multiple GPUs and show that 92% of the computation time can fully utilize the additional resources.

Waves in the Earth's core. II. Magneto–Coriolis modes

We investigate the normal modes of a rapidly rotating, electrically conducting, inviscid fluid sphere in the presence of a background magnetic field. Neglecting the fast modes that are essentially slightly modified inertial waves, we focus on the slower branch of modes traditionally called magneto–Coriolis (MC) modes. We focus on the magnetostrophic limit, an approximation that neglects the fluid inertia for these modes by setting the magnetic Rossby number E η equal to zero, converting the momentum equation from a prognostic to a diagnostic equation. In contrast to a local analysis and contrary to previous assertions, this inertialess-inviscid double limit is perfectly well behaved and shows that E η → 0 and E η = 0 are essentially identical. With applied axisymmetric background fields, we find the modes are critically damped, with decay rates greater than or equal to the eigenfrequency, in contrast to the local theory. We find evidence for both eastward and westward propagating columnar modes, in contradiction of commonly perceived wisdom. The results bode well for numerical time integrations based on the magnetostrophic approximation, which is highly suited to the core of the Earth and holds the potential for more realistic solutions of the dynamo equations.

Hydroelastic lumps in shallow water

Hydroelastic solitary waves propagating on the surface of a three-dimensional ideal fluid through the deformation of an elastic sheet are studied. The problem is investigated based on a Benney–Luke-type equation derived via an explicit non-local formulation of the classic water wave problem. The normal form analysis is carried out for the newly developed equation, which results in the Benney–Roskes–Davey–Stewartson (BRDS) system governing the coupled evolution of the envelope of a carrier wave and the wave-induced mean flow. Numerical results show three types of free solitary waves in the Benney–Luke-type equation all of which are predicted by the BRDS system: plane solitary wave, lump (i.e., fully localized traveling waves in three dimensions), and transversally periodic solitary wave. They are linked together by a dimension-breaking bifurcation where plane solitary waves and lumps can be viewed as two limiting cases, and transversally periodic solitary waves serve as intermediate states. The stability and interaction of solitary waves are investigated via a numerical time integration of the Benney–Luke-type equation. For a localized load moving on the elastic sheet with a constant speed, it is found that there exists a transcritical regime of forcing speed for which there are no steady solutions. Instead, periodic shedding of lumps can be observed if the forcing moves at speed in this range.

Dynamic calibration of differential equations using machine learning, with application to turbulence models

We present a methodology for calibration of parametric ordinary and partial differential equation models, using off-the-shelf software for back-propagation in Neural Networks (NN). As a prototypical example, we consider calibration of a Reynolds-averaged Navier-Stokes (RANS) turbulence closure model, against ground truth data from direct numerical simulations (DNS) of two different turbulent flows. Numerical time integration is represented as a custom NN, where only the RANS model parameters are trainable. A loss function is defined to quantify the mismatch between the NN prediction and the ground truth over a predefined, finite time integration window. This loss function is then minimized using a gradient descent method utilizing the back-propagation algorithm. This dynamic approach to training is to be contrasted with a static approach, wherein a least square regression estimate for parameters is obtained in the limit of an infinitesimal time integration window. In a first test of static and dynamic approaches against ground truth data generated by the model, the former proves to be significantly faster and more accurate than the latter at recovering the parameters. When both calibration approaches are tested against DNS data, for which it is known that the model cannot achieve a perfect fit, the static approach yields a good prediction only for short times, while the dynamic approach results in physical and stable predictions over the entire integration window. After optimization of the dynamic approach for time step, spatial resolution, stability, and physics-based constraints, we obtain a 50% improvement of outcomes over those obtained from the existing, manually calibrated set of parameters, demonstrating the merits of this systematic and automated procedure.

Numerical integration algorithms and constraint formulations for an ALE-ANCF cable element

In this study, the Absolute Nodal Coordinate Formulation (ANCF) is examined in the framework of the Arbitrary Lagrangian–Eulerian (ALE) formulation for cable simulation. The objective is to determine the stability, efficiency, and accuracy in different frequency domains for two different numerical time integration methods. The methods being compared are the Runge–Kutta method RK4, and the Generalized-α method. For simulations where the lower frequency domain is of interest, the Generalized-α method is stable, accurate and can be used with large time-steps; while where the high frequency domains is of interest, the RK4 procedure can be a more efficient method.The penalty method, the Lagrange multiplier method and linear coupling are tested as constraint methods for connecting ALE-ANCF cable elements. For RK4, the penalty constraint method highly affects the maximum time-step, which again affects the total solution time. For the Generalized-α method, the constraint methods all provide accurate results, and they have similar solution time. ALE-ANCF cable elements can be connected by different constraint techniques without impeding the ability to exhibit large deformations.

Large Deformation Hermitian Finite Element Coupled Thermoelasticity Analysis of Wave Propagation and Reflection in a Finite Domain

In the present paper, a finite element nonlinear coupled thermoelasticity formulation is presented for analysis of the wave propagation, reflection, and mixing phenomena in the finite length isotropic solids. The governing equations are derived based on the second Piola-Kirchhoff stress and the full form of Green’s strain-displacement tensors to account for the large deformations and finite strains. In contrast to the available researches, the assumption of very small temperature changes compared to the reference temperature is released in the present research. Galerkin’s method, a weak formulation, and cubic elements are employed to obtain the time-dependent non-linear finite element governing equations. The proposed solution procedure to the resulting highly nonlinear and time-dependent governing equations employs an updating algorithm and Newmark’s numerical time integration method. The wave propagation and reflection phenomena are investigated for both the mechanical and thermal shocks and time variations of distributions of the resulting displacements, temperature rises, and stresses are illustrated graphically and discussed comprehensively. Furthermore, the effects of the non-linear terms are discussed comprehensively. Results reveal that in the non-linear analysis, no fixed speed of wave propagation can be defined.

An efficient algorithm to compute the exponential of skew-Hermitian matrices for the time integration of the Schrödinger equation

We present a practical algorithm to approximate the exponential of skew-Hermitian matrices up to round-off error based on an efficient computation of Chebyshev polynomials of matrices and the corresponding error analysis. It is based on Chebyshev polynomials of degrees 2, 4, 8, 12 and 18 which are computed with only 1, 2, 3, 4 and 5 matrix–matrix products, respectively. For problems of the form exp(−iA), with A a real and symmetric matrix, an improved version is presented that computes the sine and cosine of A with a reduced computational cost. The theoretical analysis, supported by numerical experiments, indicates that the new methods are more efficient than schemes based on rational Padé approximants and Taylor polynomials for all tolerances and time interval lengths. The new procedure is particularly recommended to be used in conjunction with exponential integrators for the numerical time integration of the Schrödinger equation.

On a hybrid approximation concept for self‐excited periodic oscillations of large‐scale dynamical systems

AbstractWhen approximating periodic solutions in the context of large‐scale dynamical systems involving strong local nonlinearities, efficiency is of special interest. Hence, the literature suggests a combination of two approximation methods for increasing the ratio of computational cost to accuracy. Within this contribution, a combination of Finite Difference and Harmonic Balance method is proposed. Due to the usage of Harmonic Balance it is shown, that the resulting equations only depend on the degrees of freedom that are affected by nonlinear forces. As an application, a self‐excited limit cycle of a chain of oscillators is approximated and results are compared against numerical time integration to highlight qualitative accuracy.

BSTAB: an open-source software for computing the basin stability of multi-stable dynamical systems

The pervasiveness of multi-stability in nonlinear dynamical systems calls for novel concepts of stability and a consistent quantification of long-term behavior. The basin stability is a global stability metric that builds on estimating the basin of attraction volumes by Monte Carlo sampling. The computation involves extensive numerical time integrations, attractor characterization, and clustering of trajectories. We introduce bSTAB, an open-source software project that aims at enabling researchers to efficiently compute the basin stability of their dynamical systems with minimal efforts and in a highly automated manner. The source code, available at https://github.com/TUHH-DYN/bSTAB/, is available for the programming language Matlab featuring parallelization for distributed computing, automated sensitivity and bifurcation analysis as well as plotting functionalities. We illustrate the versatility and robustness of bSTAB for four canonical dynamical systems from several fields of nonlinear dynamics featuring periodic and chaotic dynamics, complicated multi-stability, non-smooth dynamics, and fractal basins of attraction. The bSTAB projects aims at fostering interdisciplinary scientific collaborations in the field of nonlinear dynamics and is driven by the interaction and contribution of the community to the software package.

Modified G-L thermoelasticity theory for nonlinear longitudinal wave in a porous thermoelastic medium

This paper presents a nonlinear analysis of the transient response of a porous medium affected by external traction. The formulation is presented based on a new modified G-L theory of thermoelasticity while temperature and strain rate parameters are included in the governing equations. Based on the finite strain theory (FST), the Lagrangian strain–displacement equation and Second Piola-Kirchhoff stress are applied to express the generalized form of thermoelasticity equations including large deformation. Both mechanical and thermal shocks are applied and finally, the effects of loading rate on transient response are discussed. An updated finite element method and Newmark’s numerical time integration method are employed to solve the nonlinear and time-dependent equations. It is determined that the modified G-L model is more powerful in predicting the wave propagation phenomenon. Furthermore, the findings indicate that the external surface shock induces compressive stresses, temperature raise, and volume fraction field increase in a porous medium. It is found compared to the classic model, the modified G-L model is superior in capturing the perceived wave propagation phenomenon.

Guaranteed globally injective 3D deformation processing

We extend recent advances in the numerical time-integration of contacting elastodynamics [Li et al. 2020] to build a new framework, called Injective Deformation Processing (IDP), for the robust solution of a wide range of mesh deformation problems requiring injectivity. IDP solves challenging 3D (and 2D) geometry processing and animation tasks on meshes, via artificial time integration, with guarantees of both non-inversion and non-overlap. To our knowledge IDP is the first framework for 3D deformation processing that can efficiently guarantee globally injective deformation without geometric locking. We demonstrate its application on a diverse set of problems and show its significant improvement over state-of-the-art for globally injective 3D deformation.

Three-dimensional aperiodic biped walking including the double support phase using LIPM and LPM

This paper presents a trajectory generation algorithm with which a three-dimensional (3D) biped robot can perform aperiodic gaits by modifying only a small set of gait parameters. In addition to the double support phase (DSP), the gait can transit smoothly from one single support phase (SSP) to another. We decouple the sagittal and coronal dynamics, firstly. In the sagittal plane, the linear inverted pendulum model (LIPM) is used to generate the walking reference trajectory during the SSP. An extra template model, linear pendulum model (LPM), is added to describe the motion in the DSP. For the coronal plane, we only utilize the LPM for trajectory generation. Thanks to linearity properties, the proposed method can obtain the biped locomotion computationally fast without the need for numerical time-integration. Herein, we also introduce the trajectory generation algorithm for different aperiodic gaits including from standing to walking, stopping walking, as well as speed switch. A full-dynamics 3D humanoid robot is used for the tests of the developed reference trajectory through simulation. The results are promising for implementations.

Explicit finite element method for nonlinear flutter analysis of composite panels

This work presents an efficient explicit finite element model for predicting the nonlinear aeroelastic behavior of composite panels in the supersonic regime. The first-order shear deformation plate theory in conjunction with the von Kármán nonlinear strains is used for structural modeling and the linear piston theory is used to model the aerodynamic loads. In order to reduce the computational cost of the simulations, a lumping procedure is employed in the mass and aerodynamic damping matrices of the finite element model. No modal reduction is performed and the central difference method is used for the numerical direct integration in time of the nonlinear equations. The model is verified using results from the literature and it is demonstrated that the lumping procedure drastically reduces the computational cost of the simulations.

Development of new embrittlement trend curve based on Japanese surveillance and atom probe tomography data

• A new embrittlement trend curve based on Japanese surveillance data was developed. • A simple analytical solution of the embrittlement trend curve (ETC) can be obtained from the original differential equation. • The new ETC can predict microscopic and macroscopic changes caused by neutron irradiation. We have developed a new embrittlement trend curve (ETC) for neutron-irradiated reactor pressure vessel(RPV) steels based on Japanese surveillance and material test reactor data. Our model is based on experimental data obtained by atom probe tomography (APT) and Charpy V-notch impact tests. The square root of the volume fraction ( V f ) of the solute atom cluster (SC) formed by neutron irradiation, which is determined by APT, is correlated to the Charpy transition temperature shift ( ΔT 41J ). In addition to the ΔT 41J and V f data, other APT data (the average volume of solute clusters (SC), the number density of SC, and the matrix Cu content in RPV steels) were also used to calibrate the coefficients in the new ETC model. The new ETC model has a simple functional form and, unlike the current version of JEAC4201, does not require numerical time integration. The new ETC model has similar ΔT 41J prediction performance to the current JEAC4201 model. Furthermore, the APT data are also predicted well, which suggests a good physical basis for the new ETC.

Nonlinear energy sink to enhance the landing gear shimmy performance

The shimmy phenomenon is a significant concern in aircraft landing gear dynamics. The prediction of the shimmy instability is an essential issue in landing gear design to develop a passive or active suppression method. This work investigates the application of the nonlinear energy sink (NES) concept to mitigate the effects of shimmy in landing gears. The NES concept has been used in recent research on mechanical vibrations. It comprises a passive target energy transfer method that refers to a one-way energy transfer from a primary to a nonlinear subsystem. The landing gear model is based on torsional displacement coupled with the tyre classical elastic string analogy model. The NES device connects to the wheel shaft, and it comprises a mass, a linear damper, and a pure cubic spring. The numerical integration in time was used to assess the shimmy onset speed and the post-shimmy limit cycle oscillations. A parametric analysis of the landing gear nonlinear dynamics without the NES is presented. The design space of possible NES parameters is given, obeying design constraints, and inspected to assess adequate NES designs. The best NES samples are included in the landing gear dynamics to study their influence. Results have shown that the NES can adequately expand the operational speed range for no shimmy and lead to lower LCO amplitudes in the post-shimmy for a reasonable range of speeds. The NES concept’s successful employment for the landing gear dynamics suggests an enormous potential form of passive shimmy control.

Investigation of critical delamination characteristics in composite plates combining cubic spline piezo-layerwise mechanics and time domain spectral finite elements

The scope of the present work is to elucidate the reverberations of propagating guided waves through delaminated regions confined in composite plate structures. As it is well known, delamination between consequent plies is the most common type of damage in composite structures. A mean of investigation of delamination cracks is through guided wave propagation. It has been observed that either type of guided waves (antisymmetric or symmetric) has different repercussions inside and outside the delaminated area. Inside this work is used a well-established numerical tool developed by the authors. By making good use of its reported numerical capabilities, the authors used the Hermitian Cubic Splined Layerwise mechanics to incorporate the delamination and the physical presence of piezoelectric actuators. The two mentioned aspects were then integrated into a time domain spectral finite element. The outcome was the development of a multinode explicit layerwise finite element, able of modelling the physical presence of active piezoelectric sensors, the existence of delamination cracks through the thickness and the excitation of each fundamental wave mode separately (A0 &S0). The calculated semi-diagonal mass matrix, due to the collocation of the nodes with the grid points of Gauss–Lobatto–Legendre quadrature, boosts the speed of numerical time integration, while the use of high order polynomials as shape functions, conduce to a vast reduction of the mesh size, enabling the thorough investigation and simulation of high frequency wave propagation.