Efficient and robust implementation of the TLISMNI method

Abstract The dynamics of large scale and complex multibody systems (MBS) that include flexible bodies and contact/impact pairs is governed by stiff equations. Because explicit integration methods can be inefficient and often fail in the case of stiff problems, the use of implicit numerical integration methods is recommended in this case. This paper presents a new and efficient implementation of the two-loop implicit sparse matrix numerical integration (TLISMNI) method proposed for the solution of constrained rigid and flexible MBS differential and algebraic equations. The TLISMNI method has desirable features that include avoiding numerical differentiation of the forces, allowing for an efficient sparse matrix implementation, and ensuring that the kinematic constraint equations are satisfied at the position, velocity and acceleration levels. In this method, a sparse Lagrangian augmented form of the equations of motion that ensures that the constraints are satisfied at the acceleration level is used to solve for all the accelerations and Lagrange multipliers. The generalized coordinate partitioning or recursive methods can be used to satisfy the constraint equations at the position and velocity levels. In order to improve the efficiency and robustness of the TLISMNI method, the simple iteration and the Jacobian-Free Newton−Krylov approaches are used in this investigation. The new implementation is tested using several low order formulas that include Hilber–Hughes–Taylor (HHT), L-stable Park, A-stable Trapezoidal, and A-stable BDF methods. The HHT method allows for including numerical damping. Discussion on which method is more appropriate to use for a certain application is provided. The paper also discusses TLISMNI implementation issues including the step size selection, the convergence criteria, the error control, and the effect of the numerical damping. The use of the computer algorithm described in this paper is demonstrated by solving complex rigid and flexible tracked vehicle models, railroad vehicle models, and very stiff structure problems. The results, obtained using these low order formulas, are compared with the results obtained using the explicit Adams–Bashforth predictor-corrector method. Using the TLISMNI method, which does not require numerical differentiation of the forces and allows for an efficient sparse matrix implementation, for solving complex and stiff structure problems leads to significant computational cost saving as demonstrated in this paper. In some problems, it was found that the new TLISMNI implementation is 35 times faster than the explicit Adams–Bashforth method.


Introduction
The numerical solution of constrained MBS problems requires the numerical integration of differential and algebraic equation (DAE) systems. Potra [1] and Negrut [2] discussed several techniques for the numerical solution of the DAE system in MBS dynamics. Large scale and complex MBS applications that include flexible bodies and contact/impact elements lead to stiff problems. Because explicit integration methods can be very inefficient and often fail in the case of stiff problems, the use of implicit numerical integration methods is recommended. A good understanding of the process of converting the DAE system to a second order ordinary differential equations (ODE) system, along with the ability to efficiently generate the needed derivative information for the implicit integration, allowed developing robust implicit numerical methods [3,4]. Nonetheless, existing implicit numerical integration methods used for the solution of MBS applications have several drawbacks. First, most of these methods do not ensure that the nonlinear algebraic constraint equations are satisfied at all levels (position, velocity, and acceleration). Second, existing implicit methods require the use of numerical differentiation of the force vectors in order to solve a nonlinear system of algebraic equations using a Newton-Raphson algorithm; analytical differentiation cannot in general be used with general MBS algorithms that allow for the use of tabulated and Spline data to define the forces and constraints.
Third, many of the existing implicit MBS integration methods are not suited for an efficient sparse matrix implementation because they lead to dense matrices. In order to address these limitations, the two-loop implicit sparse matrix numerical integration (TLISMNI) method was proposed [5]. The TLISMNI method ensures that the constraint equations are satisfied at the position, velocity, and acceleration levels, does not require the use of numerical or analytical differentiation of the forces, and allows for efficient sparse matrix implementation.
In the case of constrained MBS dynamics, the algebraic kinematic constraint equations are used to describe mechanical joints and specified motion trajectories that restrict relative and absolute motion of the MBS components. Different techniques have been proposed in the literature for the numerical integration of DAE systems. These techniques include the direct integration method, the generalized coordinates partitioning method, and the constraint stabilization method [3,6]. The direct integration method employs a numerical integration method of ordinary differential equations to integrate the differential algebraic equations without any modification of the integration algorithm or dynamic equations. This integration method, which integrates the second time derivative of the generalized coordinates without considering their dependency, is simple, easy to implement, and computationally fast, but it suffers from a lack of error control on the constraints and may lead to erroneous results [7]. In the generalized coordinates partitioning method, proposed by Wehage and Haug [8], the system generalized coordinates are partitioned into dependent and independent coordinates, and only the independent accelerations are integrated to define the state of the system at the next time step. Dependent coordinates are obtained by solving the nonlinear algebraic constraint equations using a Newton-Raphson algorithm. This method has the advantage of ensuring that the constraint equations are satisfied to within a user specified tolerance. A good prediction of the dependent generalized coordinates can improve the efficiency of the method, since such a prediction can satisfy the constraint equations without the need for using Newton-Raphson iterations.
Numerical experimentation has shown that, with a good prediction of the dependent variables, Newton-Raphson iterations are needed only for a few time steps during the dynamic simulation.
Baumgarte [9] proposed a different approach based on a constraint stabilization method. In this approach, the constraint equations and their derivatives are used with the second derivative of the constraint equations to form a penalty term. Using this penalty term, the constraint stabilization integration algorithm requires direct integration of all the accelerations. The disadvantage of this method, however, is that there is no general and uniformly accepted method for selecting the coefficients of the constraint equations and their derivatives in the penalty expression. An improper selection of these coefficients can lead to wrong solutions and make the method less robust.
As previously mentioned, the implicit integration methods can be more efficient than explicit integration methods in solving stiff DAE systems [10]. Implicit methods require the solution of a system of nonlinear algebraic equations at each time step. Newton's method or modified Newton's method appears to be the most widely used approach for iteratively solving the resulting implicit method nonlinear algebraic equations. For large problems, most of the calculations required for the implicit integration are associated with the computation of the Jacobian and the solution of the resulting large system of nonlinear equations. This is in addition to the need for a relatively large core memory for the storage of the coefficient matrix and its decomposition. Gear and Saad [11] proposed the use of Krylov-subspace projection method known as the incomplete orthogonalization method (IOM) and Arnoldi's algorithm, which are iterative methods for the solution of linear systems [12]. Arnlodi's algorithm and IOM do not require the factorization of the coefficient matrix in any form, and therefore, less storage as compared to the direct methods is needed. Brown and Hindmarsh [13] referred to the combined stiff ODE method and Krylov method as a matrix-free method and discussed the theoretical and computational aspects of the combined algorithm. Brown and Hindmarsh [13] viewed the combined Newton-IOM and Newton-Arnoldi method as an inexact Newton method, which belongs to a class of methods in which the linear system of Newton iteration is solved approximately. The implementations of the Arnoldi's algorithm and IOM for solving linear system of equations are discussed in detail by Saad [12]. This paper discusses TLISMNI implementation issues, including the step size selection, the error control, and the effect of the numerical damping and proposes a new efficient and robust TLISMNI-based solution procedure for constrained MBS equations. Furthermore, the paper examines the use of the generalized coordinate partitioning and recursive methods in the TLISMNI framework solution; both approaches are used in order to satisfy the constraint equations at the position and velocity levels. In both cases, the constraints are also satisfied at the acceleration level by using the augmented Lagrangian form to solve for the accelerations and Lagrange multipliers. The use of the computer implementation described in this paper is demonstrated by solving complex rigid and flexible body tracked vehicle models, railroad vehicle models, and very stiff structural problems. Several second order formulas such as Hilber-Hughes-Taylor (HHT) method, including numerical damping, L-stable Park formula, Astable Trapezoidal formula, and A-stable BDF formula are used in this investigation as the integration formulas, and recommendations are made with regard to the appropriateness of each of these integration formula for a particular MBS application. The efficient integration of the Krylov subspace projection method with these integration formulas within the TLISMNI framework solution procedure is one of the main contributions of this investigation. The results, obtained using these integration methods, are compared with the results obtained using the explicit Adams predictor-corrector method. The TLISMNI method does not require numerical differentiation of the forces, allows for an efficient sparse matrix implementation for solving complex and very stiff structure problems, and ensures that the constraint equations are satisfied at the position, velocity, and acceleration levels. The use of this method significantly improves the computational efficiency as demonstrated in this paper using several examples. In some stiff problems, it was found that the new TLISMNI implementation is 35 times faster than the explicit Adams-Bashforth method.

Background
In this section, the constrained MBS equations of motion and the solution procedures used to solve these equations are briefly discussed. The definitions provided in this section will be used repeatedly in this paper.

MBS equations of motion
In MBS dynamics, the constraint relationships are used with the differential equations of motion to solve for the unknown accelerations and constraint forces. While this approach leads to a sparse matrix structure, it has the drawback of increasing the problem dimensionality and requiring more sophisticated numerical algorithms to solve the resulting DAE system. Using the generalized absolute Cartesian coordinates, the equations of motion of a body i can be written as [7,14] where i M is the mass matrix of the body, The symbols that appear in this equation ( As previously mentioned, i k F and i k T are the generalized joint forces associated, respectively, with the translation and orientation coordinates of body i . Using the results of Eq. (3), the reaction forces at the joint definition points can be determined using the concept of the equipollent systems of forces.

Generalized coordinate partitioning
In order to ensure that the algebraic kinematic constraint equations are satisfied at the position and velocity levels, the independent accelerations i q  are identified and integrated forward in time in order to determine the independent velocities i q  and independent coordinates i q .
Knowing the independent coordinates from the numerical integration, the dependent coordinates d q can be determined from the nonlinear constraint equations using an iterative Newton-Raphson algorithm that requires the solution of the system  converge. On the other hand, the generalized coordinate partitioning technique is suited for the sparse matrix implementation and can be used for MBS applications that include rigid and flexible bodies and systems that suffer from singularity problems, such as closed chains.

Algorithm and sparse matrix implementation
In order to ensure that the constraint equations are satisfied at the position level, the dependent coordinates are determined by solving the nonlinear algebraic constraint equations using the iterative Newton-Raphson procedure. In the sparse matrix implementation one can use the following system of algebraic equations in Newton iterations [7,8]: In this equation, q C is the constraint Jacobian matrix, q is the vector of Newton differences, and d I is a Boolean matrix that has zeroes and ones only; with the ones in the locations corresponding to the independent coordinates in order to ensure that i  Δq 0. The square coefficient matrix in Eq. (4) is sparse, and therefore, sparse matrix techniques can be used to efficiently solve the preceding system of equations for the dependent coordinates. Once the dependent coordinates are determined, the dependent velocities can be obtained using the following linear sparse system that defines the constraint equations at the velocity level: The right hand side of Eq. (5) is assumed to be known since i q  is determined from the numerical integration, and t C depends on time and the coordinates that are assumed to be known from the position analysis. One advantage of the structure of Eqs. 4 and 5 is that if the set of independent coordinates change during the simulation time one has only to change the locations of the nonzero entries of the matrix d I , while the structure of the Jacobian matrix q C remains the same. Once the generalized coordinates and velocities are determined, the augmented form of Eq.
(2) can be constructed and solved for the acceleration q  , and Lagrange multipliers λ as previously mentioned. The main steps for a numerical algorithm using the generalized coordinate partitioning can then be summarized as follows: 1. An estimate of the initial conditions that define the MBS initial configuration is made.
The initial conditions that represent the initial coordinate and velocities must be a good approximation of the initial configuration of the system.
2. Using the coordinates, the constraint Jacobian matrix q C can be constructed, and an LU factorization algorithm can be used to identify a set of independent coordinates.
3. Using the values of the independent coordinates, the constraint equations   ,t  C q 0 can be considered as a nonlinear system of algebraic equations in the dependent coordinates.
This system can be solved iteratively using Eq. (4) and a Newton-Raphson algorithm that employs sparse matrix techniques. 7. If the end of the simulation is not reached, the algorithm returns to Step 2.
In the following section, it is shown how the steps of this generalized coordinate partitioning solution procedure can be modified for the TLISMNI implementation discussed in this paper.

TLISMNI algorithm
The TLISMNI method was recently proposed for the solution of the MBS differential/algebraic equations [15,16]. As mentioned before, this method ensures that the algebraic constraint equations are satisfied at the position, velocity, and acceleration levels; does not require numerical or analytical differentiation of the forces, and allows for efficient sparse matrix implementation. In this section, the solution steps for a modified TLISMNI algorithm that allows for the implementation of different low order integration formulas are presented. The performance of this algorithm will be tested in a later section of the paper using complex flexible and rigid body vehicle models. The two-loop implicit procedure proposed in this study can be algorithm that employs the generalized coordinate partitioning can be summarized as follows: 1. Assuming that the state of the system is known at time n t , the constraint Jacobian matrix q C can be constructed, and used with Gaussian elimination procedure to determine a set of system independent coordinates i q . Therefore, the vector of system generalized coordinates q can be partitioned to dependent and independent coordinates as where d q is the set of system dependent coordinates. , and ( 1) , 1 7. Using the convergence and error criteria introduced in a later section of this paper, one can judge whether or not the integration is successful. If the solution satisfies the convergence and error criteria, the coordinates, velocities, accelerations, and Lagrange multipliers are accepted. In this case, one needs to calculate the new step size in order to advance the integration. If, on the other hand, convergence is not achieved or the error exceeds the specified error tolerance, then more iterations are required using a time step size that achieves convergence. In this case, the algorithm returns to Step 2 and the process continues until the convergence and error criteria are satisfied.
8. This process continues until the desired end of the simulation time is reached.
It is clear from the computational algorithms presented in this section that the TLISMNI method uses sparse coefficient matrices that require minimum storage. Furthermore, no numerical differentiation of the external or inertia forces with respect to the coordinates and velocities is required. The algorithm outlined above is equivalent to performing fixed point iteration for the solution of nonlinear system of equations for 1 n q , and 1 n q  using an implicit integration formula [17,18]. The TLISMNI algorithm employs the implicit integration formulas which can be more efficient than the explicit formulas in many applications, particularly in the case of stiff problems.

Krylov subspace and inexact Newton method
One of the main contributions of this paper is to implement the Krylov subspace projection method in the TLISMNI algorithm. In the TLISMNI algorithm described in this paper, one can consider the outer loop as a set of ODE initial value problem that can be written as . Assuming these ODE's are stiff and there is a need to use an implicit integration formula to obtain a solution, one can write the general formula for implicit integration methods as where ψ is a vector that has variables previously computed, h is the time step size,  is a constant that depends on the implicit formula used, and f y is the right hand side of the system of differential equations. In this paper, an equivalent form of Eq. (6) will be used. This form can be written in terms of Typically, Newton and modified Newton algorithms are employed to solve Eq. (7) If we consider  J 0, the iterative procedure will be simple and is referred to in this case as simple iterations, and it will be equivalent to performing fixed point iteration for the solution of the nonlinear equations for 1 n y as follows: Equation (9) represents the outer loop for the algorithm described in Section 3. In case  J 0 expensive computations are required for the numerical calculations or approximations of J .
Furthermore, matrix decomposition, and large storage space will be required, especially when solving large scale problems. In addition, the numerical differentiation required for a general MBS implementation can be a source of numerical errors that may lead to non-accurate solutions as the problem dimensionality increases. Therefore, for such complex and large scale problems, the use of a method that can approximately solve Eq. (8) at reasonable cost and also reduces the core memory required is desirable. The Incomplete orthogonalization method (IOM) and Arnoldi's algorithm are examples of such methods that will be discussed in a following subsection. IOM and Arnoldi's algorithm are methods for the approximate solution of a linear system Ax = b in N R [12]. The use of such methods in the solution of the nonlinear system represents the amount by which the vector ( ) m s fails to satisfy the Newton equation (Eq. (8)). The theoretical foundation and convergence of the inexact Newton method is discussed by Brown and Hindmarsh [13].

Krylov subspace projection method
In this subsection, the iterative Arnoldi's algorithm for the solution of linear system Ax = b is briefly reviewed [12]. For the nonlinear problem given in Eq. (7), which can be written as κ . More details about the method can be found in [12,13]. The Arnoldi's algorithm can be described as follows [12]: All the features and properties of Alrnoldi's algorithm still hold for IOM [12]. ( 1 m  ). Therefore the approximate solution can be given as Due to the structure of m U , the vector m p can be calculated using the previous i p 's and m v as follows: Similarly, because of the structure of m L , the vector m z can be updated using An extensive literature about the Jacobian-free Newton-Krylov method can be found in [19]. It can be shown for the algorithms described in the preceding subsection that the matrix A is not needed explicitly. One needs to compute the matrix-vector product Av only. Since for the stiff ODE, it is assumed that ( ) h     A F x I J , where x is an approximation to the root of ( )  F x 0 , then the matrix-vector product Av in the above algorithm can be replaced by different quotients of the form where  is a scalar. The resulting algorithm can be referred to as a finite-difference projection method or Jacobian-free Newton Krylov Method.

TLISMNI Newton-Krylov algorithm implementation
In this section, the use of Jacobian-free Newton-Krylov approach for solving stiff DAE is presented. The proposed algorithm will employ the scaled Arnoldi's or IOM method [13], where the weight associated with component i y of y during the iteration is where RTOL and ATOL are the relative and absolute error tolerances, respectively. In order to 1. Assuming that the state of the system is known at time n t , the constraint Jacobian matrix q C can be constructed, and used with Gaussian elimination procedure to determine the set of system independent coordinates i q . Therefore, the vector of system generalized coordinates q can be partitioned to dependent and independent coordinates as where d q is the set of dependent coordinates. 7. If the convergence criteria proposed in the following section is satisfied and the error is less than the user specified tolerance, the coordinates, velocities, and the solution for the acceleration and Lagrange multipliers are accepted. In this case, one needs to update the history, and calculate the new step size in order to advance the integration. If convergence is not achieved or the error exceeds the specified error tolerance, then the time step should be reduced and the algorithm is restarted. In this case, the algorithm goes to step 2 again until the convergence and error criteria are satisfied.
8. This process continues until the desired end of the simulation time is reached.
The proposed algorithm takes advantage of the Jacobian-Free Newton-Krylov algorithm, and still exploits the sparse matrix structure to achieve minimum storage. Furthermore, the constraints are satisfied at the position, velocity, and acceleration levels, and no numerical differentiation of the external or inertia forces with respect to the coordinates and velocities is required to obtain the Jacobian matrix.

Convergence criteria
The TLISMNI method can be designed to have an iterative outer loop to solve for the coordinates, velocities, accelerations, and Lagrange multipliers. The steps of the TLISMNI outer loop can be summarized as follows: 1-Define the system coordinates 0 The algorithm outlined above is equivalent to performing fixed point iteration for the solution of nonlinear equations for 1 n q , and 1 n q  using an implicit integration formula [17,18] More details on the linear convergence analysis can be found in [13]. The condition used in this section implies that the iteration converges in a sufficient small ball around the root. It can be shown that TLISMNI's iteration is much less expensive than the Newton iterations and allows for much more rapid variation of the step size in order to achieve convergence. In order to ensure rapid convergence in a practical implementation, the convergence rate C is selected to be much smaller than 1, for example 0.3. The goal is to have an algorithm that takes less than four iterations to converge, particularly in the case of large scale problems. On the other hand, this restriction can slow TLISMNI simple iteration method and the implicit formula can lose its properties, such as numerical damping in case of HHT [15]. In such a case, it is recommended to use the Jacobian-Free Newton-Krylov method [18] instead of the TLISMNI simple iteration method. The convergence procedure for both the TLISMNI simple iteration method and the TLISMNI Newton-Krylov method can be described as [13]:

Low order integration formulas
As previously mentioned, implicit integration methods transform the MBS differential equations of motion into a set of nonlinear algebraic equations. These nonlinear algebraic equations can be solved iteratively for the required solution. Explicit methods can be very inefficient or fail to solve stiff problems which are characterized by widely separated eigenvalues because of high frequency contents. Higher order integration formulas cannot be used effectively for the solution of such stiff problems. For these problems, low order A-stable integration formulas can be more efficient. Using A-stable low order integration formula, the restriction on the size of the time step to maintain absolute stability is no longer required. In this section, different low order integration formulas will be discussed and recommendations are made on the appropriateness of each method for a particular problem.

Park method
The first integration formula considered in this section is the Park method, which has order 2 and was proposed by Park [20] as an improved stiffly stable method for direct integration of nonlinear structural dynamic equations. By combining the Gear's two-step and three-step method, a superior stiffly stable method was developed [20,21]. Park method can be applied to both stiff and non-stiff problems. The results indicate that Park method is second best after the trapezoidal rule for non-stiff problems and appears to be stable for stiff problems that include frictional contact/impact phenomena, as will be demonstrated in the numerical results section.
Pogorelov [22] proposed the use of Park method for the solution of stiff constrained MBS applications. Park method does not require any history derivative information, which can cause numerical instability in nonlinear dynamics problems even though the methods are unconditionally stable. The equations that define the generalized coordinates and velocities at time 1 n t  in Park method are given, respectively, by where h is the time step, and subscript n refers to vectors at time n t . The local truncation error of Park method in terms of the accelerations is The second integration method considered in this investigation is the second order backward differentiation formula (BDF2) method. This method was proposed by Gear [23]. The BDF2 method is a stiffly A-stable method that has been widely used for the solution of stiff problems due to their good stability properties for such problems. The BDF2 equations for the generalized coordinates and velocities are and 1 1 1 As is clear from these equations, the BDF2 method does not require any history derivative information; therefore, it will be stable with stiff problems as is the case with Park method [20].
The local truncation error for the BDF2 method is Comparing the BDF2 method's truncation error and the Park method's truncation error, one can see that Park method can achieve the same accuracy as the BDF2 method for a larger time step size by a factor of approximately 1.5.

Trapezoidal method
The third integration method that will be considered in this investigation is the trapezoidal method, which is considered the most accurate second order A-stable methods. The method does not damp out any frequency content from the system and it is unconditionally stable for linear problems, and conditionally stable for nonlinear systems [17,24]. More details about the stability and the accuracy of the trapezoidal method can be found in [25]. The trapezoidal method equations do require history derivative information; therefore the method suffers from stability problems with stiff problems. The trapezoidal algebraic equations that represent the generalized coordinates and velocities are given, respectively, as The truncation error using the trapezoidal method can be estimated as where  and  are constants. The differential equations of motion and the constraint equations can be written in the form The Nemark method is a first order method that produces a second order accurate method when The HHT method introduces a numerical damping parameter  to the equations of motion in order to allow for energy dissipation and retain the order and stability condition of the method.
The modified equations of motion can be written as In order to obtain a stable solution using the implicit HHT method, the following relations should be satisfied 0.3 0 The local truncation error using HHT method is

Error control and time step selection
The error criteria and the time step-size selection for the proposed TLISMNI algorithm are discussed in this section. Using the estimate of the truncation errors where s is a safety factor and h is the current time step size. In case the time step size does not satisfy the user specified error, e   , then the time step size is reduced and the iteration step is restarted.

Explicit Adams method
Adams method is an explicit predictor-corrector method for the numerical solution of first order ordinary differential equations. This method cannot be used to directly solve a system of differential and algebraic equations. The algebraic equations must be first eliminated using the generalized coordinate partitioning technique, which is in principle equivalent to the embedding technique [7,14,27]. The embedding technique, that eliminates the reaction forces and the dependent coordinates, leads to a minimum set of differential equations expressed in terms of the independent coordinates (degrees of freedom) only. Using the generalized coordinate partitioning, one can write the system coordinates as where i q is the vector of independent coordinates or degrees of freedom, and d q is the vector of dependent coordinates.
One can rewrite the equations of motions in terms of the independent coordinates as follows [7,14]: In this equation, is the velocity transformation matrix, and Note that the augmented formulation previously discussed in this paper can be used to determine all the accelerations. The independent accelerations can then be identified and integrated forward in time. Therefore, the explicit Adams method can be used with both the augmented formulation and the embedding technique.
In order to use the explicit Adams method, one has to transform the second order ordinary differential equations to a system of first order ordinary differential equations (ODE). This is accomplished by introducing the state vector Taking the time derivative of this vector and substituting the value of i q  from the solution of Eq. (30) yields The right hand side of this equation is a function of i q , i q  andt . Therefore, the original problem is reformulated as   ,t  y f y  , which is the standard form of a first order ODE that can be solved by Adams method. The explicit Adams method used in this study is the Adams predictor-corrector method documented in the book by Shampine and Gordon [26]. In this method, a predicted value ,t f y , that is, One can then use the corrected value 1 n y to evaluate the function 1 n f which is used in the next time step. A full documentation of the procedures used in this explicit method for error control and selection of the order and time step can be found in the literature [27].

Numerical examples
In this section, numerical results obtained using a simple pendulum, a tracked vehicle model, and railroad vehicle models, are used to demonstrate the use of the TLISMNI algorithm proposed in this investigation for solving large and complex stiff systems that include flexible bodies and contact/impact forces. The results obtained using the TLISMNI algorithm and the explicit Adams predictor-corrector method, are compared in terms of efficiency and accuracy. In addition, recommendations are made on the appropriateness of each integration formula for a particular problem.

Pendulum example
The pendulum used in this example is assumed to be initially horizontal and fall under the effect of the gravity forces. The pendulum model is developed using the finite element absolute nodal coordinate formulation (ANCF). The beam in this pendulum is discretized using twodimensional finite beam elements along its length as shown in Fig. 1 [28]. The pendulum is assumed to have undeformed length 0.4 m, cross sectional area 0.04 0.04  2 m , the mass density 7200 kg/m 3 , modulus of elasticity 11 2 2 10 N/m  , and Poison's ratio 0.3. In order to compare the performance of the TLISMNI algorithm and the explicit Adams method, a dynamic simulation of the stiff flexible pendulum is performed by using the two integration methods. Figure 2 shows that the implicit TLISMNI method damped out some high frequency oscillations and this integration method produces a smoother solution than the one obtained using the explicit Adams method. With regard to the efficiency, the TLISMNI method is 35 times faster than the Adams method for this example. Also it is important to mention that there is no significant difference between the two reference motion solutions, and therefore, the numerical damping of the TLISMNI method does not have a significant effect on the rigid body motion of the beam, as will be demonstrated using more complex examples.

Rigid tracked vehicle model
The tracked vehicle model shown in Fig. 3 is a challenging model that represents an armored personnel carrier consisting of a chassis, 2 idlers, 2 sprockets, 10 road-wheels, and 128 track links (64 for each track). Figure 5 shows the engagement of the track links with some of the vehicle components. The vehicle has a suspension system that consists of road arms placed between the road wheels and the chassis as well as shock absorbers connected to each road arm, as shown in Fig. 4. Table 1 shows the stiffness and the damping coefficients of the contact models and the suspension system used for this model. The road arms and the sprockets are  be at least six times faster than the explicit Adams method on the same machine. It is also important to mention that the recursive approach was found to be more efficient than the augmented formulation for such a complex chain problem.

Flexible tracked vehicle model
In this model, the vehicle model shown in Fig. 3 is used, where a new compliant continuumbased joint formulation is used for the joint formulation between the track links. In the numerical investigation presented in this paper, a three-dimensional cable element is used to model the flexibility of the chain links [29,30].

Summary and conclusions
The objective of this paper is to integrate the Newton- TLISMNI implementation issues including step size selection, convergence criteria, error control, and the effect of the numerical damping were discussed. Simple pendulum, complex rigid and flexible tracked vehicle, and railroad vehicle models were used to demonstrate the use of the proposed TLISMNI implementation. A comparison between the results obtained using the TLISMNI algorithm and the explicit Adams predictor-corrector method showed good agreement.
On the other hand, using TLISMNI method, which does not require numerical differentiation of the forces and allows for an efficient sparse matrix implementation for solving complex and very stiff structure problems, significantly improves the simulation time. For the rigid body model considered in this investigation, the TLISMNI is at least five times faster than the explicit Adams method. Using the TLISMNI algorithm with integration formulas that employ numerical damping such as HHT in the simulation of flexible body models considered in this study can achieve up to thirty five times faster simulation compared to Adams method. Nonetheless, it is important to mention that there are cases of non-stiff problems in which the use of explicit Adams method can be more efficient than the TLISMNI methods. The use of the Jacobian-Free Newton-Krylov approach instead of the simple iteration approach improves the convergence and accuracy of the TLSMNI method. Preconditioning and parallelization techniques need to be explored in a computational framework based on the TLISMNI Newton-Krylov approach, which improves the convergence for the stiff equations. Also the use of the generalized coordinate partitioning proved to be efficient in the case of differential/algebraic equation problems. More investigations are required in order to develop better criteria for the selection of the best set of independent coordinates and the time to change these coordinates. This is necessary in order to improve the performance of the TLISMNI method and avoid singularity problems arising in applications that include closed chains.