Application analysis of time stepping methods for impact problems

Dynamic contact and impact problems are widely applicable. An accurate solution method for these kinds of problems can be used in many fields of mechanical engineering (e.g., cutting metalwork, cogwheel drives, etc.). However, the proper handling of the contact is problematic, as there emerges a substantial amount of nonlinearity in the displacement field. Therefore, a spurious high frequency oscillation is present in the solution. These oscillations must be avoided, as divergence can easily occur in the contact algorithm due to them. In order to eliminate this effect, the applied numerical method must be chosen and set properly. In this study, a comprehensive guide is provided for the appropriate selection of the proper numerical method and its parameters.

and the Hilber-Hughes-Taylor-α (HHT-α) method [16] will be studied, which are widely used backward increment schemes. The main advantage of these methods is that they are unconditionally stable. This is a major aspect; hence these methods are often built in commercial FE software. The Newmark method is the most common implicit method in this kind of software. The formulation of this method contains a parameter through which the numerical damping of the method can be affected directly. The HHT-α method is the improved version of the Newmark method, where the dissipation level can be controlled with two parameters. Nevertheless, the solution produced with a backward increment scheme is more time-consuming to achieve, as a system of equations must be solved in every time step. There are many more time integration schemes, like the Houbolt, the Wilson-θ, the generalized-α, etc. that can be found in [9]. Besides, there exist several multi stage schemes, like Kim's method [17] and the Kolay-Ricles-α (KR-α) method [18]. In this paper, only the above mentioned four methods are considered, as they proved to be the most useful for contact problems, based on our experiences.
The right choice of the applied numerical method is crucial, as it has a significant impact on the solution characteristics, the computation time and the achievable accuracy. There is not any universal recipe to reach an ideal solution, many different schemes must be taken into consideration to examine, which method gives the best solution as of the user's preferences. It is a further problem that due to the spatial discretization a spurious high frequency oscillation occurs in the resulting fields of displacement and contact pressure, which distorts the correct solution. This difficulty can be at least partially eliminated by applying a certain amount of numerical damping by the solution of the equation of motion. Nevertheless, it is also substantial to filter only the disturbance out by leaving the useful component of the signal intact. Thus, the adjustability of the numerical damping is an important aspect to consider, when choosing the right method. The main difficulty is to find an optimal region for the dissipation parameter to achieve a solution with acceptable accuracy, while significantly reducing the spurious oscillations.
In this article, the effectiveness of the aforementioned kinds of methods is going to be analyzed. The accuracy and the practical applicability of the examined methods will be studied from several aspects. Besides, parameter identifications are going to be made in order to determine the right amount of numerical damping for methods with adjustable dissipation level. The main goal of this paper is to provide a comprehensive application analysis about how to choose the right method with the correct parameters for a dynamic contact problem.

TIME STEPPING METHODS
Let us consider a structural dynamic problem, where the interacting bodies are linearly elastic. The semi-discretized equation of motion after spatial discretization by finite element method can be written as a system of linear differential equations of second order according to where M, C and K denote the mass, damping and stiffness matrices, respectively. On the right side f t means the applied load vector, while u t , _ u t and € u t are the unknown displacement, velocity and acceleration vectors at time instant t, respectively. Moreover the solution of this differential system of equations requires both initial and boundary conditions from which the former can be written as where u 0 and v 0 denote the initial displacement and velocity vectors, respectively. In the following, four different numerical methods will be introduced for solving numerically the equation of motion (1), hence the deduction will be presented for each method in the next four chapters.

Central difference method
The CDM is a commonly used forward increment method for various mechanical applications. The main benefits of using this scheme are that it is second order accurate, and computations can be executed very fast due to its explicit formulation. However, the central difference method is only conditionally stable, which reduces the applicable time step region. In this method, the velocity and the acceleration are given using the displacement at the current, previous and next time step the following way.
By substituting Eqs (3) and (4) into the Eq. of motion (1) and reordering it, Eq. (5) can be obtained. The main problem with this formula is that it requires the inversion of the term ΔtC þ 2M, which is a computationally expensive procedure, Hence, it is worthwhile to combine Eq. (4) with Eq. (6), as it results a formulation that only requires the inversion of M. As the inversion of a lumped mass matrix is computationally cheap, it is rewarding to diagonalize the mass matrix by the method described in [12]. It is important to note that the formula applied for the _ u t velocities originates from the backward Euler method. Thus, the applied explicit technique is a certain kind of combination of the central difference and the backward Euler method, Equations (4) and (6) can be substituted into the equation of motion (1). By ordering it the same way, as it was evaluated for Eqs (5) and (7) can be obtained.

Backward Euler method
The backward Euler scheme has an implicit character, so, a system of equations (which is often nonlinear) needs to be solved for every time step, resulting high computational costs. This method is unconditionally stable, first order accurate and has a certain amount of numerical damping. The main problem is, that the dissipation level is not controllable; hence it often leads to energy dissipation of an undesirable rate, The general formulation of this method according to [13] can be written as Eqs (8) and (9). Eq. (1) of motion for time step t þ Δt leads to formulation (10).

Newmark method
The Newmark method [14,15] is a commonly used implicit integration scheme, which is built in many commercial FEM software (Ansys, Abaqus, etc.), Its basic formulation is determined by Eqs (11) and (12). Based on [14], the optimal parameter setting can be given as Eq. (13). By setting the δ parameter according to Eq. (13), the integration accuracy and stability can be obtained by only one parameter. By changing this parameter, the amount of numerical damping also alters; hence the problematic oscillations are fully or partially extinguishable. This method is unconditionally stable and first order accurate (except by γ 5 0.5 which results second order accuracy). The main disadvantage is the lack of high-frequency damping which limits the reduction of spurious noises.

HHT-α method
The HHT-α method [16] is another scheme, which is widely employed for solving the equation of motion of mechanical systems. This method is the improved version of the Newmark method with controllable numerical dissipation through the α parameter. The main benefit in relation to the Newmark scheme is that by this method numerical damping can be exerted without degrading the order of accuracy, The basic Equations (11) and (12) of the Newmark method remain valid for this scheme. However, the equation of motion (10) must be modified using the α parameter according to Eq. (14). In the interval of α ∈ [À0.5, 0], the amount of numerical damping can be controlled where the smaller value of α means the smaller amount of dissipation, If α is treated independently, then Eq. (13) parameters from the Newmark method remain unchanged, so the result of the numerical computation can be influenced by two separate parameters, δ and α. However, if β and γ are chosen according to Eq. (15), the scheme becomes second order accurate, although this does not cause the improvement of the numerical results necessarily. In this case the practical range α ∈ ½−1=3 ; 0. In practice, it is usually worthwhile to try both Eq. (15) and the independent choice of δ and α .

FEM model
The methods that have been described above should be tested using an example by which it can be judged, how well they can handle the sudden changes in characteristic indicators like the velocity or the contact pressure. For this purpose, a 1D impact problem is an ideal choice as it is a widely used test example in the literature of contact problems [11,17,18]. The main benefits of this example are that its solution can be made easily and fast and it has an exact analytical solution, which provides a good reference to assess the results. Besides, it challenges the achievable accuracy of the numerical method at a sufficient level. Thus, if an accurate solution can be produced for this test example, then the applied method can presumably treat more complex geometries too. Therefore, a linearly elastic rod is considered which is moving towards a rigid wall (Fig. 1).
In order to solve this test example, a FEM model must be constructed based on the mechanical model shown in Fig. 1. This FEM model is built using 1D truss elements, which possess 2 nodes with 1 Degree of Freedom (DoF) per node, which is the longitudinal displacement. The FEM model is shown in Fig. 2, where L ¼ 20 ½mm; h ¼ 0:1 ½mm; n ¼ 200 ½ − ; and v 0 ¼ 1000 ½mm=s. The model contains uniform elements with E ¼ 90 ½MPa and ρ ¼ 7:85$10 −9 ½kg=mm 3 material parameters. Based on the performed calculations, above the applied element number the results do not improve substantially.

Treatment of the contact
The proper handling of the contact is an important aspect by the modeling.
As it has been mentioned earlier, by the Lagrange multiplier method, the penetration between the contacted bodies is completely prohibited. Because of this condition, certain displacement constraints are prescribed the following way: Gðu þ XÞ ¼ 0: (16) In the equation above, X means the initial configuration (vector of length n þ 1 filled with the initial coordinates of nodes), u is the displacement vector, and G is the contact constraint matrix. This equation is a certain kind of manifestation of the Hertz-Signorini-Moreau (HSM) conditions [4] in contact mechanics. The detailed description of here applied formulation can be found in [11]. Thus, Eq. (1) must be modified with regard to the contact force acting between the elastic rod and the rigid wall So, the equation of motion can be written in the form of Eq. (17) where λ t denotes the Lagrange multiplier at time instant t.
In this case, that quantity is equivalent with the surface contact pressure. Thus, the Lagrange multiplier method proceeds at the here presented example by treating λ as unknown and solving Eqs (16) and (17) simultaneously. The numerical methods examined in Section 2 can easily be deducted for the contacting case by replacing the equation of motion (1) with (17).

PARAMETER IDENTIFICATION
In the case of the Newmark and the HHT-α method, the quality of the numerical solution highly depends on the proper choice of the alterable parameter(s). Hence, a parameter identification is needed in order to test the applicability for a contact problem, Thus, it is practical to define certain indicators by which it can be quantified, how accurate the 1D example may be solved and how efficient the elimination of the spurious oscillations is. Firstly, relative errors are considered, which can be calculated by Eq.

Newmark method
Based on these two indicators the practical range of δ parameter for a contact problem may be ascertained. After determining them for the contact pressure for multiple values of δ inside the interval δ ∈ ½0 ; 1, these indicators can be plotted in the function of δ according to Fig. 3a and b, respectively. As it can be seen in these figures, NoP p is decreasing strictly monotonic, while ε p is starting to grow after a certain value of δ is stepped over. Hence, an ideal parameter value cannot be ascertained, by which both of the applied indicators are reaching their minimum, so an acceptable compromise must be made. Based on Fig. 3, a so called "practical stage" may be assigned, inside, which both measures remain at sufficiently low values. In accordance with the numerical tests using the FEM model, this practical stage can be determined as α ∈ ½0:3 ; 1.

HHT-α method
The parameter identification for this method is more complicated, as there are two independent parameters affecting the numerical solution. The applied indicators are the same as by the Newmark method, but by reason of the one more ruling parameter, a surface plot must be applied ( Fig. 4a and b). In order to establish a beneficial correlation between α and δ, the Period Error (PE) [19] function must be analyzed, Based on this examination, Eq. (19) may be concluded to reach a parameter setting that results minimum for the period error. Considering Fig. 4, this connection proves to be a practical compromise between the accuracy and the elimination of the spurious oscillations. Hence, it is beneficial to choose the parameter values according to Eq. (19), especially as it makes the parameter setting much easier. Based on Fig. 4 and in accordance with the concrete numerical tests with the FEM model, the "practical stage" can be ascertained as α ∈ ½−0:5 ; − 0:3.

APPLICATION IN A 1D EXAMPLE
In order to compare the applied methods and check the results of the parameter identification, a concrete function must be plotted that characterizes the 1D impact. As the spurious oscillations are mainly problematic in the contact pressure and the velocity functions, the methods shall be compared using the time evolution of these quantities (Figs 5-6). The magnified area on the left side shows the oscillations emerging directly after the moment of impact. The enlarged region on the right side of Fig. 5 indicates how steep the pressure decreases to zero which illustrates the accuracy of the solutions. The prominence in the zoomed plot on the right side of Fig. 6 is caused by the wave propagation in the rod and must be as narrow as possible.
Besides, the indicators that have been used in section 4 are illustrated by identical dissipation for each method in Figs 7 and 8 for the pressure and the velocity, respectively. By studying them, it can be assessed that there cannot be found an ideal method which to use is beneficial in all aspects.
The central difference method offers a fast and accurate result, but the large amount of oscillations makes it problematic to apply. In contrast, the backward Euler method offers the best performance in the elimination of oscillations, but this feature is coupled with the highest relative error. Based on this study, the Newmark and the HHT-α prove to have the highest potential at this kind of application. Both of them are controllable with one or more parameters, by which an acceptable compromise may be made. Even from these two methods, the HHT-α seems to possess more beneficial damping characteristics, as by identical dissipation it results a better solution. As the Newmark method is built in the ANSYS mechanical APDL software, the performed calculations were verified by that comparison. Applying the same FEM model with identical parameters in ANSYS,   the results show a very precise match, proving that the calculations made in Julia language are accurate.

CONCLUSION
In this study a thorough analysis was made about the application of numerical methods in contact problems. Four different, widely used numerical schemes were deducted and applied for a 1D example. In the parameter identification, different aspects were presented, through which a correct parameter setting may be determined for numerical methods with one or more control parameters. The minimization of the period error proved to be a beneficial compromise by the HHT-α method regarding the examined measures. The results of the application in a concrete 1D example also confirmed that this approach is correct. Nevertheless, a numerical method that is ideal to solve a contact problem in all aspects cannot be ascertained, as the right choice highly depends on the user's actual priorities. The main function of this paper is to provide a comprehensive guide for the correct selection of the numerical method for a concrete contact problem.