Implicit Time Integration Method Research Articles

A multirate time integrator for regularized Stokeslets

The method of regularized Stokeslets is a numerical approach to approximating solutions of fluid–structure interaction problems in the Stokes regime. Regularized Stokeslets are fundamental solutions to the Stokes equations with a regularized point-force term that are used to represent forces generated by a rigid or elastic object interacting with the fluid. Due to the linearity of the Stokes equations, the velocity at any point in the fluid can be computed by summing the contributions of regularized Stokeslets, and the time evolution of positions can be computed using standard methods for ordinary differential equations. Rigid or elastic objects in the flow are usually treated as immersed boundaries represented by a collection of regularized Stokeslets coupled together by virtual springs which determine the forces exerted by the boundary in the fluid. For problems with boundaries modeled by springs with large spring constants, the resulting ordinary differential equations become stiff, and hence the time step for explicit time integration methods is severely constrained. Unfortunately, the use of standard implicit time integration methods for the method of regularized Stokeslets requires the solution of dense nonlinear systems of equations for many relevant problems. Here, an alternate strategy using an explicit multirate time integration scheme based on spectral deferred corrections is incorporated that in many cases can significantly decrease the computational cost of the method. The multirate methods are higher-order methods that treat different portions of the ODE explicitly with different time steps depending on the stiffness of each component. Numerical examples on two nontrivial three-dimensional problems demonstrate the increased efficiency of the multi-explicit approach with no significant increase in numerical error.

Fast simulation code for heating, phase changes and dopant diffusion in silicon laser processing using the alternating direction explicit (ADE) method

Enhancements to our existing finite-differences code for the simulation of laser heating, melting and evaporation of silicon are presented. The Knudsen evaporation model has been added to the previously used enthalpy-based model in order to simulate laser pulses with pulse lengths down to a few nanoseconds. Fick’s diffusion law has also been incorporated allowing laser doping by dopant diffusion in silicon melt to be described. Finally, the basic equations for the alternating direction explicit method (ADE) have been adapted to consider nonlinear temperature-enthalpy relations, thus including affects due to phase changes. This improved the simulation speed by up to factor of 100 compared to standard explicit and implicit time integration methods. Details of the ADE algorithm and numerical stability issues are presented in this paper. Validation of the code is presented by comparing to different integration methods and to experimental results. The final code successfully simulates melting, evaporation and dopant diffusion by multiple laser pulses in three dimensions in an acceptable computing time.

Assessment of ecosystem health of tropical shallow waterbodies in eastern India using turbulence model

In the present study, a numerical model of the hydrodynamic and thermal structure of artificial shallow lakes in eastern India has been developed as a tool to assess the ecological water quality, driven by the meteorological forcings. It allows quantification of the vertical mixing processes that govern not only the thermal structure but also nutrient exchanges and the distribution of dissolved and particulate matter among water layers. Vertical temperature profiles were calculated by solving coupled partial differential equations for temperature (heat energy balance equation), one-dimensional momentum equation and a second order closure scheme for small-scale turbulence effects, i.e. turbulent kinetic energy and turbulent dissipation rates numerically using an implicit time integration method. The effect of advection due to the inflow and outflow is not taken into consideration as these shallow waterbodies are assumed to behave as a closed lake. The oxygen level measured at different layers during the simulation period in these waterbodies reflects hypolimnetic oxygen depletion due to the thermal stratification in the aquatic environment. The changes in the stratification regime in these waterbodies are expected to affect the water quality and the health of the ecosystem, primarily based on temperature and dissolved oxygen parameter and in particular, the certain features of the oxygen resources of the hypolimnion.

Viscous nonlinearity in central difference and Newmark integration schemes

Of the various implicit and explicit time integration methods the central difference method and Newmark’s method are the most popular ones to solve Cauchy’s equation of motion in structural solid dynamics. While in their standard forms elastic and viscous linearity are presumed, they are applicable to nonlinear elastic materials with little alterations. However, the case of viscous nonlinearity is rarely investigated, although it demands some algorithmic modifications of the methods in concern. In the context of a finite element spatial discretization we will provide a description of those emended integration schemes and illustrate them exemplarily.

Efficient Fully Implicit Time Integration Methods for Modeling Cardiac Dynamics

Implicit methods are well known to have greater stability than explicit methods for stiff systems, but they often are not used in practice due to perceived computational complexity. This paper applies the backward Euler (BE) method and a second-order one-step two-stage composite backward differentiation formula (C-BDF2) for the monodomain equations arising from mathematically modeling the electrical activity of the heart. The C-BDF2 scheme is an L-stable implicit time integration method and easily implementable. It uses the simplest forward Euler and BE methods as fundamental building blocks. The nonlinear system resulting from application of the BE method for the monodomain equations is solved for the first time by a nonlinear elimination method, which eliminates local and nonsymmetric components by using a Jacobian-free Newton solver, called Newton--Krylov solver. Unlike other fully implicit methods proposed for the monodomain equations in the literature, the Jacobian of the global system after the nonlinear elimination has much smaller size, is symmetric and possibly positive definite, which can be solved efficiently by standard optimal solvers. Numerical results are presented demonstrating that the C-BDF2 scheme can yield accurate results with less CPU times than explicit methods for both a single patch and spatially extended domains.

Implementation of an ANCF beam finite element for dynamic response optimization of elastic manipulators

Theoretical and practical aspects of an absolute nodal coordinate formulation (ANCF) beam finite element implementation are considered in the context of dynamic transient response optimization of elastic manipulators. The proposed implementation is based on the introduction of new nodal degrees of freedom, which is achieved by an adequate nonlinear mapping between the original and new degrees of freedom. This approach preserves the mechanical properties of the ANCF beam, but converts it into a conventional finite element so that its nodal degrees of freedom are initially always equal to zero and never depend explicitly on the design variables. Consequently, the sensitivity analysis formulas can be derived in the usual manner, except that the introduced nonlinear mapping has to be taken into account. Moreover, the adjusted element can also be incorporated into general finite element analysis and optimization software in the conventional way. The introduced design variables are related to the cross-section of the beam, to the shape of the (possibly) skeletal structure of the manipulator and to the drive functions. The layered cross-section approach and the design element technique are utilized to parameterize the shape of individual elements and the whole structure. A family of implicit time integration methods is adopted for the response and sensitivity analysis. Based on this assumption, the corresponding sensitivity formulas are derived. Two numerical examples illustrate the performance of the proposed element implementation.

Soil-Water Coupled Finite Deformation Analysis Based on a Rate-Type Equation of Motion Incorporating the SYS Cam-Clay Model

This paper presents a new method of soil-water coupled finite deformation analysis of saturated soils that considers inertial forces. This method allows changes in the geometric shape of the soil to be taken into account and is capable of dealing with all types of external forces irrespective of whether they are static or dynamic. To be more specific, the paper describes the following points, which differ from the conventional methods: 1) the governing equations for saturated soil including the rate-type equation of motion containing a jerk term of the soil skeleton conforming to u-p formulation and updated Lagrangian, 2) derivation of a weak form of the rate-type equation of motion and discretization of the finite elements, and 3) use of the implicit time integration method for application of the conventional linear acceleration method (which assumes linear variation of acceleration) to the jerk term. By mounting the elasto-plastic constitutive equation (SYS Cam-clay model), which can cover a wide range of soils and soil conditions, onto the above method of analysis, examples of simulation of dynamic/static triaxial laboratory testing of saturated soil specimens are described. The soil specimens were assumed to be medium dense sand under conditions of small-amplitude cyclic loading, partial drainage, and constant cell pressure. The simulation yielded the following results: (1) In the case of low frequencies, compaction occurs during loading and compression progresses over the entire specimen. (2) In the case of high frequencies, during loading and in the period in which wave propagation continues within the specimen after the end of loading, compaction occurs at the drained end of the specimen, whereas liquefaction occurs in its interior. After this stage, massive compression takes place within the specimen, leading to consolidation (consolidation after liquefaction).

Nonlinear dynamic analysis of an eccentrically prestressed damped beam under a concentrated moving harmonic load

This paper focuses on the geometrically nonlinear dynamic analysis of an eccentrically prestressed simply supported damped beam subjected to a concentrated moving harmonic load. The nonlinear dynamic deflections of the beam are obtained by polynomial functions. The Kelvin–Voigt model for the material of the beam is used. Two coupled nonlinear systems of equations of motion are derived using Lagrange's equations under the assumptions of the Euler–Bernoulli beam theory with the von-Kármán's nonlinear strain–displacement relationships. The rotary inertia, axial displacement and axial inertia are included in the formulation. The nonlinear equations of motion are solved by using the implicit time integration method of Newmark- β in conjunction with the Newton–Raphson iteration method. In this study, the effects of large deflections, the internal damping of the beam, the velocity of the moving harmonic load, the prestress load, the eccentricity of the prestress load and the excitation frequency on the dynamic response of the beam are discussed. The obtained results are compared with the results based on the linear beam theory. Convergence studies are performed. Numerical results show that the above-mentioned effects play a very important role in the deflections of the beam.

Evaluation of effective parameters in metal bellows forming process

Metal bellows have wide applications in aerospace, micro-electromechanical and industrial systems. Forming process of the metal bellows is very sensitive to increasing the ratio of crown to root diameter. In this state, precise control of the parameters is very important in order to form high-quality metal bellows with good thickness distribution and desirable dimensions and resilience. In this paper, a new method has been proposed for manufacturing of the metal bellows and important parameters such as initial length of tube, internal pressure, axial feeding and velocity, mechanical properties and the type of materials were investigated by finite element (FE) analysis (LS-Dyna) and experimental tests. The explicit time integration method is used for modeling the tube-bulging and folding processes. Meanwhile, the implicit time integration method is used for the spring back stage. Finally, the results of finite element method (FEM) and experiments show a very good agreement. The results of the present work could be used as a basis of designing a new type of the metal bellows.

Continuum Dual-scale Modeling of Liquid Composite Molding Processes

In this study, a comprehensive continuum dual-scale model has been developed for LCM processes for woven fiber reinforcements. This model incorporates the main flow mechanisms involved, including mass exchange between dual-scale porosities, micro-diffusion within fiber tows between neighboring repeating unit cells (RUCs), and the capillary effect within fiber tows. Macro-and micro-flows are simulated at the same numerical scale. A sink and a source term are introduced into the two different scaled flows to couple them together. An unsaturated flow model is adopted to describe micro-diffusion within fiber tows between neighboring RUCs. The physics of the sink/ source term is modelled according to the impregnation of tows within a RUC. Two models developed in soil physics have been adopted to evaluate the capillary effect and the partially saturated permeability within fiber tows. A finite element method is employed to solve the governing equations, and a simple implicit time integration method is suggested. A 1D RTM example has been modeled. The modeling results indicate that all of the three main mechanisms have significant influences on the LCM processes, but the mass exchange between the dual porosities plays a dominate role, particularly under high injection rates.

Non-linear vibrations of laminated cylindrical shallow shells under thermomechanical loading

The geometrically non-linear vibrations of linear elastic composite laminated shallow shells under the simultaneous action of thermal fields and mechanical excitations are analysed. For this purpose, a model based on a very efficient p-version first-order shear deformation finite element, with hierarchical basis functions, is employed. The equations of motion are solved in the time domain by a Newmark implicit time integration method. The model and code developed are partially validated by comparison with published data. Parametric studies are carried out in order to study the influence of temperature change, initial curvature, panel thickness and fibre orientation on the shells’ dynamics.

Unconditionally Stable Fully Explicit Finite Difference Solution of Solidification Problems

An unconditionally stable fully explicit finite difference method for solution of conduction dominated phase-change problems is presented. This method is based on an asymmetric stable finite difference scheme for approximation of diffusion terms and application of the temperature recovery method as a phase-change modeling method. The computational cost of the presented method is the same as the explicit method per time-step, while it is free from time-step limitation due to stability criteria. It robustly handles isothermal and nonisothermal phase-change problems and is very efficient when the global temperature field is desirable (not accurate front position). The robustness, stability, accuracy, and efficiency of the presented method are demonstrated with several benchmarks. Comparison with some stable implicit time integration methods is also included. Numerical experiments show that the evaluated temperature field for large values of the Fourier numbers has good/reasonable agreement with the result of small Fourier numbers and the exact/reference solution.

On the robustness of a multigrid method for anisotropic reaction-diffusion problems

In this paper, we consider a reaction-diffusion boundary value problem in a three-dimensional thin domain. The very different length scales in the geometry result in an anisotropy effect. Our study is motivated by a parabolic heat conduction problem in a thin foil leading to such anisotropic reaction-diffusion problems in each time step of an implicit time integration method [7]. The reaction-diffusion problem contains two important parameters, namely e >0 which parameterizes the thickness of the domain and μ >0 denoting the measure for the size of the reaction term relative to that of the diffusion term. In this paper we analyze the convergence of a multigrid method with a robust (line) smoother. Both, for the W- and the V-cycle method we derive contraction number bounds smaller than one uniform with respect to the mesh size and the parameters e and μ.

Numerical analysis of time integration errors for nonequilibrium radiation diffusion

Numerical analysis of time integration errors for nonequilibrium radiation diffusion is considered. Two first-order implicit time integration methods are studied. Asymptotic analysis and modified equation analysis are applied to both time integration methods. Numerical experiments are used to highlight the results of the analysis. Asymptotic analysis is used to highlight the source of temperature spiking when a hot radiation wave propagates into a cold material. Modified equation analysis is used to provide insight into the thermal wave speed coming from the two different first-order methods.

Linearly implicit time integration methods in real-time applications: DAEs and stiff ODEs

The methods for the dynamical simulation of multi-body systems in real-time applications have to guarantee that the time integration of the equations of motion is always successfully completed within an a priori fixed sampling time interval, typically in the range of 1.0–10.0 ms. Model structure, model complexity and numerical solution methods have to be adapted to the needs of real-time simulation. Standard solvers for stiff and for constrained mechanical systems are implicit and cannot be used straightforwardly in real-time applications because of their iterative strategies to solve the nonlinear corrector equations and because of adaptive strategies for stepsize and order selection. As an alternative, we consider in the present paper noniterative fixed stepsize time integration methods for stiff ordinary differential equations (ODEs) resulting from tree-structured multi-body system models and for differential algebraic equations (DAEs) that result from multi-body system models with loop-closing constraints.

Unconditionally Stable Fully Explicit Finite Difference Solution of Solidification Problems

An unconditionally stable fully explicit finite difference method for solution of conduction dominated phase-change problems is presented. This method is based on an asymmetric stable finite difference scheme for approximation of diffusion terms and application of the temperature recovery method as a phase-change modeling method. The computational cost of the presented method is the same as the explicit method per time-step, while it is free from time-step limitation due to stability criteria. It robustly handles isothermal and nonisothermal phase-change problems and is very efficient when the global temperature field is desirable (not accurate front position). The robustness, stability, accuracy, and efficiency of the presented method are demonstrated with several benchmarks. Comparison with some stable implicit time integration methods is also included. Numerical experiments show that the evaluated temperature field for large values of the Fourier numbers has good/reasonable agreement with the result of small Fourier numbers and the exact/reference solution.

Stability of semidiscrete formulations for elastodynamics at small time steps

Solutions of direct time-integration schemes for elastodynamics that converge in time to conventional semidiscrete formulations may be polluted at small time steps by spurious oscillations that violate the principle of causality, for example by arising before wave fronts. This degradation is not an artifact of the time-marching scheme, but rather a property of the solution of the semidiscrete formulation itself. An analogy to singularly perturbed elliptic problems provides an upper bound on the time step for the onset of these oscillations. A simple procedure of spatial stabilization is proposed to remove this pathology from implicit time-integration schemes, without affecting unconditional temporal stability. Spatially stabilized implicit time-integration methods are free of noncausal oscillations at small time steps.

HLI, A Direct Method Suitable for Partial and Fully Implicit Time Integration of Primitive Equation Meteorological Models

HLI, the homogeneously linearized implicit method is a direct implicit time integration method for nonlinear prognostic equations. For homogeneous linear test problems it is identical to the known implicit method, when applied to the full model area. As most theoretical results concerning the implicit method use linear homogeneous test problems, these apply immediately to HLI. The implicit method, when applied to the whole model area, is called the fully implicit method. In order to reduce the numerical expense for large integration areas, a partial implicit method is defined by posing artificial numerical boundary conditions, which result in a smaller set of equations, corresponding to the part of the integration area chosen. The application of either full or partial methods to nonlinear and inhomogeneous problems is achieved by defining for each grid point the local homogeneously linearized problem. Partial schemes are convenient when using multiprocessor computers, as only local and neighbouring data have to be used at each grid point. Only grid point discretisations of finite order are admitted. The proposed method would not be suitable for the spectral method. Computational examples are presented to show the suitability of this method for inhomogeneous situations and the ease of the treatment of internal boundary conditions. A semi-implicit treatment was also tested, which treats only the terms responsible for the fast waves implicit. This requires only the Fourier transformation of one field and in the generalised Fourier back transformation only one linear equation has to be solved, rather that a system of equations. In this paper only the trapezoidal implicit method was used. HLI could, however, be applied in the same way with other implicit methods.

Thermally induced transitions to chaos in plate vibrations

The geometrically nonlinear vibrations of linear elastic and isotropic plates under the combined effect of thermal fields and mechanical excitations are analysed. With this purpose, a model based on a p-version, hierarchical, first-order shear deformation finite element is employed. The equations of motion are solved in the time domain by Newmark's implicit time integration method. The temperature and the amplitude of the mechanical excitation are varied, and transitions from periodic to non-periodic motions are found.

Analysis of Thin Beams and Cables Using the Absolute Nodal Co-ordinate Formulation

The purpose of this paper is to present formulations for beam elements based on the absolute nodal co-ordinate formulation that can be effectively and efficiently used in the case of thin structural applications. The numerically stiff behaviour resulting from shear terms in existing absolute nodal co-ordinate formulation beam elements that employ the continuum mechanics approach to formulate the elastic forces and the resulting locking phenomenon make these elements less attractive for slender stiff structures. In this investigation, additional shape functions are introduced for an existing spatial absolute nodal co-ordinate formulation beam element in order to obtain higher accuracy when the continuum mechanics approach is used to formulate the elastic forces. For thin structures where bending stiffness can be important in some applications, a lower order cable element is introduced and the performance of this cable element is evaluated by comparing it with existing formulations using several examples. Cables that experience low tension or catenary systems where bending stiffness has an effect on the wave propagation are examples in which the low order cable element can be used. The cable element, which does not have torsional stiffness, can be effectively used in many problems such as in the formulation of the sliding joints in applications such as the spatial pantograph/catenary systems. The numerical study presented in this paper shows that the use of existing implicit time integration methods enables the simulation of multibody systems with a moderate number of thin and stiff finite elements in reasonable CPU time.