A generalized approach for implicit time integration of piecewise linear/nonlinear systems

A generalized solution scheme using implicit time integrators for piecewise linear and nonlinear systems is developed. The piecewise linear characteristic has been well‐discussed in previous studies, in which the original problem has been transformed into linear complementarity problems (LCPs) and then solved via the Lemke algorithm for each time step. The proposed scheme, instead, uses the projection function to describe the discontinuity in the dynamics equations, and solves for each step the nonlinear equations obtained from the implicit integrator by the semismooth Newton iteration. Compared with the LCP‐based scheme, the new scheme offers a more general choice by allowing other nonlinearities in the governing equations. To assess its performances, several illustrative examples are solved. The numerical solutions demonstrate that the new scheme can not only predict satisfactory results for piecewise nonlinear systems, but also exhibits substantial efficiency advantages over the LCP‐based scheme when applied to piecewise linear systems.

others, 18,19 can be designed to realize unconditional stability for linear problems, so their time step sizes only depend on the desired solution accuracy, at the cost of greater computational effort. When applied to linear systems, implicit methods require the factorization of the effective stiffness matrix; for nonlinear systems, an iterative solution approach is inevitable. A more comprehensive review of time integration methods can be found in Tamma et al. 20 Considering piecewise linear and nonlinear problems, a simple way for numerical solutions is to employ an explicit integrator directly, like in Refs. 4,21. These methods update the current step by known state variables, and then make the updated variables satisfy the piecewise characteristics. However, as mentioned above, one drawback of explicit methods is their conditional stability, and since the dynamic equations cannot be satisfied implicitly at each time step, their solutions for nonlinear problems may be not so reliable.
On the other hand, there are also some attempts to provide available schemes using implicit integrators. Their main difficulty is to solve the implicit dynamics equations with piecewise characteristics at each time step. Yu 22 and Fadaee et al. 23 proposed that piecewise linear equations can be transformed into equivalent linear complementarity problems (LCPs), 24 and then solved by the Lemke algorithm. They use the Bozzak-Newmark method 25 as the implicit integrator in the schemes. A similar idea was employed by He et al., 26 based on the precise integration method (PIM) 27  has been used to deal with the discontinuities caused by impacts and frictions, is introduced to model the piecewise linear and nonlinear characteristics. Without loss of generality, the generalized-α method 14 is employed to integrate the dynamics equations, and the resulting nonlinear equations are solved at each time step by the semismooth Newton iteration. 32 Compared with the LCP-based schemes, the proposed approach can be directly used to solve general piecewise nonlinear systems. At the same time, the proposed approach shows significant efficiency advantages for systems with a large number of piecewise features, due to the fast convergence rate of the Newton iteration. The proposed approach is applied to several numerical examples, to assess its numerical performances. This paper is organized as follows. In Section 2, the dynamic equations for piecewise linear and nonlinear systems are presented using the projection function. The computational procedures are provided in Section 3, based on the generalized-α method and the semismooth Newton iteration. Numerical experiments are implemented in Section 4. Conclusions are finally drawn in Section 5.

| FORMULATION
For illustrative purposes, the simple piecewise linear SDOF system shown in Figure 1 is discussed first. Its elastic restoring force, N x ( ), can be formulated as follows: where x is the displacement of the cart, Δ is the initial clearance, and k 1 and k 2 are the stiffness parameters of the system. Using this simple model, the basic ideas of the LCP-based scheme are reviewed here.
From Equation (1), the complementarity relationship can be established by introducing two auxiliary variables, g and y, as  y g = 0. (2d) Using the variable y, Equation (1) can be transformed into After applying the desired time integration scheme, the numerical displacement x k at time t k satisfies k x k y r * + = * k k k 2 (4) where k* is the effective stiffness value, which is a constant for linear systems; r* k is the effective load, which needs to be updated for each time step. By combining Equations (2) and (4), one obtains the LCP y k k k g k k r k Δ = ( * + ) * + ( * + ) ( * − * ), At each step, y k and g k are first obtained by solving the LCP via the Lemke algorithm. Then x k is computed from Equation (2a) as In the proposed scheme, the projection function is introduced to formulate more general forms of piecewise linear and nonlinear characteristics. From convex analysis, the projection of a point w in a convex set C is the closest point in C to w, defined as follows: By using the projection function, Equation (1) can be transformed It indicates that the projection function in Equation (8) can replace the role of the auxiliary variable y of the LCP-based scheme in Equation (3).
For general multistage piecewise linear systems, as follows: where the constants z z z , , …, n 1 2 are defined as required to make N x ( ) continuous at x x x , , …, n 1 2 −1 , the equivalent form using the projection functions can be expressed as follows: where the convex sets are set as  0 + uniformly. Some special cases, such as where Δ 1 and Δ 2 are the initial clearances, can be transformed into the simpler form When extended to piecewise nonlinear systems, such as the equivalent form can be written as follows: The contribution of each term in Equation (14) can be decoupled as follows:  Combined with the continuity conditions, that is , it can be easily checked that Equation (14) is the same as the original Equation (13). In a similar way, for general polynomial functions, as where p p p , ,…, n 1 2 are positive integers, the equivalent form using the projection function is , Some special cases, like (18) can also be formulated as the simpler form Since formulating nonlinearities as polynomial functions can cover most situations produced in structural dynamics, other possible forms, which can also be expressed using the projection function by selecting proper To present a general formulation of piecewise linear and nonlinear MDOF systems, an auxiliary variable   y q×1 is introduced as follows: T are projected functions with respect to x. Then, the general form of the dynamic equations can be expressed as follows: where   M s s × and   C s s × are the mass and damping matrices, assumed constant for simplicity, the over dot indicates (21) constitute the dynamics equations to be solved. The introduction of the variable y is not strictly necessary, but it makes the expressions more readable. In particular, the dynamics equations of piecewise linear systems can be written as follows: and , ,

| COMPUTATIONAL PROCEDURE
Without loss of generality, the generalized-α method is employed to illustrate the computational procedures. The standard generalized-α method in Chung and Hulbert 14 only shows firstorder accuracy for acceleration, so an improved scheme in Arnold and Brüls, 33 which holds second-order accuracy for displacement, velocity, and acceleration, is adopted. It introduces an accelerationlike variable a as follows: where the subscript k denotes the state variables of step k α , , and δ are parameters of the algorithm. Using the acceleration-like variable a, the displacement and velocity of step k + 1 are updated according to is the time step size, and β and γ are additional parameters of the algorithm. The state variables of step k + 1 satisfy the discrete dynamic equations, as

Mx
Cx For piecewise linear systems, the dynamic equations to be solved are The generalized-α method embraces several well-known From linear spectral analysis, a set of optimal parameters, controlled by the spectral radius ρ ∞ at the high-frequency limit, for the generalized-α method was given as follows: A smaller ρ ∞ denotes worse accuracy and stronger numerical damping, which can be used to filter out inaccurate high-frequency contributions. Using Equation (30), the generalized-α method becomes a single-parameter scheme with respect to ρ ∞ , and to avoid excessive loss of accuracy, ρ ≥ 0.6 ∞ is recommended for general purpose integration.

| Computational procedure for piecewise linear systems
Considering the dynamic equations of piecewise linear systems in Equations (28)-(29), only Equation (29) needs to be solved iteratively, so the computational procedure is relatively simple. Substituting Equations (24)- (25) into Equation (28) yields where Since the effective stiffness matrix S is constant, S −1 can be prepared before the step-by-step solutions. Using Equation (31), Equation (29) can be rewritten as The projection function x proj ( ) C is nondifferentiable at the boundary of the convex set C. According to the semismooth Newton method, 32 the generalized derivative of x proj ( ) C can be defined as follows: Consequently, at each step, the nonlinear equation to be solved is Its Jacobian matrix with respect to y k+1 is where The iterative scheme can be written as follows: where the superscript l denotes the iteration number. In Equation (36) For clarity, the step-by-step solution using the generalized-α method for piecewise linear systems is summarized in Table 1. Other implicit integrators can be used in a similar way. Compared with the LCP-based scheme, the proposed scheme employs the semismooth Newton iteration instead of the Lemke algorithm to solve y k+1 . Since the semismooth Newton iteration has locally super-linear convergence rate, the state of the last step is used as the prediction.
Consequently, if the current time step is in a smooth interval, the problem degenerates to a linear problem, requiring only one iteration.
Only when the switching point is passed, the proposed scheme may need to use more iterations to reach a convergent solution. The computational efficiency of the proposed scheme and the LCP-based scheme is compared based on the numerical example in Section 4.

| Computational procedure for piecewise nonlinear systems
Considering the dynamic equations of piecewise nonlinear systems of Equations (26)- (27), the nonlinear equation that need to be solved at each time step is where y px x = ( ),k k k +1 +1 +1 and ẍ k+1 can be expressed as functions of x k+1 using Equation (38). The Jacobian matrix of g x ( ) k+1 with respect to x k+1 has the form The iterative scheme is In this case, matrices K K W , , x y , and Λ need to be updated at each iteration. Until the convergence condition is satisfied, x k+1 is obtained, and other state variables can be computed using Equation (38). The step-by-step solution is presented in Table 2.
LCP-based schemes have not yet been applied to piecewise nonlinear problems. Therefore, the previously described

Initialize x x x a
,˙,¨, 0 0 0 0 , and y 0 ; 3. Select the time step size h, the algorithmic parameters ρ ∞ , the tolerance error ϵ, and the maximum number of iterations N; 4. Calculate integration constants: B. For each time step 1. Calculate the effective load r k+1 : c. Update y k+1 y y Gy g y l l = + 1, d. If l N < and   g > ϵ, go to b; If l N < and   g ≤ ϵ, go to 3; If l N ≡ , abort.

Solve for
, and a k+1 : ZHANG ET AL. The generalized-α method is known as second-order accurate for smooth systems, but previous studies 34,35 show that discontinuities in velocity or acceleration caused by impacts or friction can reduce its accuracy to first-order. For this reason, it is necessary to discover the accuracy order of this method for systems with piecewise linear characteristics. Therefore, the SDOF system employed in Worden et al. 36 to model a cracked beam is solved to assess the convergence rate. The equation of motion can be written as follows: 4. Calculate integration constants: For each time step +1 and a k+1 : 4. Update x k+1 : 5. If l N < and   g > ϵ, go to 2; If l N < and   g ≤ ϵ, go to next step; If l N ≡ , abort. To show the convergence rate, the relative errors of displacement, velocity, and acceleration versus the step size h are plotted in Figure 3. Here the relative errors within [0, 10 s] are computed according to where n is the total number of steps, x k and x t ( ) k , respectively, denote the numerical and reference solutions. It can be observed that the integrators still hold second-order accuracy for displacement, velocity, and acceleration, so the piecewise characteristic does not cause the reduction of accuracy order. This is because all state variables in the recursive schemes are continuous, and no sudden change brings additional numerical errors.
The conclusion can be naturally extended to general MDOF systems, which are equivalent to the superposition of a set of SDOF systems.
Besides, the results also show that TR (ρ = 1 ∞ ) is more accurate than the dissipative schemes, so it is employed in the following examples. Despite this, the dissipative schemes are very useful in large finite element systems and stiff problems, to filter out the high-frequency dynamics and make the solutions smoother and more stable.

| MDOF case
To compare the computational efficiency of the proposed scheme with that of LCP-based ones, a set of MDOF cases, as  Using h = 0.001 s, the transient response within [0, 1 s] is tracked, and the required CPU times for these two schemes are listed in Table 3. It follows that as the number of piecewise characteristic increases, the new method gradually exhibits a significant efficiency advantage. During the numerical experiments, the Newton iteration always reaches convergence in no more than five iterations, but the number of pivoting required by the Lemke algorithm increases rapidly as the number of piecewise characteristics increases. As shown in Murthy, 37 Lemke's algorithm needs exponential complexity in the worst case, which means that it requires up to 2 q pivoting steps to reach the final result. Therefore, the proposed scheme appears to be substantially more efficient when applied to systems with a large number of piecewise features.

| Rolling mill
As shown in Figure 5, the rolling mill system, 2 a typical piecewise nonlinear problem, is considered. The governing equation can be written as follows: Using the projection function, k x ( ) can be converted to  can be used to study more complex piecewise nonlinear systems.

| MDOF case
The MDOF model in Section 4.1 is reused here, but the elastic force is replaced by Other sets and parameters remain unchanged. The initial displacements satisfy static equilibrium.
Using 10 DOFs and h = 0.0001s, Figure 8 shows the solutions of x x , 1 5 , and x 10 . As can be seen, the numerical results agree very well with the reference ones. It demonstrates that the proposed scheme works very well also for piecewise nonlinear MDOF system. For more complex problems, the convergence of the Newton iteration greatly depends on the predicted value at each step. A small time step, an accurate predicted value, and appropriate algorithmic dissipation are all helpful to improve the convergence.

| CONCLUSIONS
This study provides a solution scheme using implicit integrators for solving the dynamic response of general piecewise linear and nonlinear systems. The proposed scheme employs the projection function to express the piecewise characteristics, instead of converting them into equivalent LCPs. Compared with the existing LCP-based scheme, the new scheme possesses two benefits. First, it can be applied to general piecewise nonlinear systems, because the discrete nonlinear equations containing the projection functions can be solved uniformly using the semismooth Newton method. Additionally, according to the numerical experiment, the proposed scheme is more efficient than the LCPbased scheme for the system with many piecewise features. The computational procedures are presented using the generalized-α method and the semismooth Newton iteration. Other implicit integrators can also be utilized in a similar manner. Numerical experiments demonstrate that the Newton iteration is more efficient than the Lemke algorithm, especially when the system has a lot of piecewise linear features. When applied to piecewise nonlinear problems, the numerical solutions of the new scheme present a high consistency with the reference ones.