Application of a Novel Picard-Type Time-Integration Technique to the Linear and Non-Linear Dynamics of Mechanical Structures: An Exemplary Study

: Applications of a novel time-integration technique to the non-linear and linear dynamics of mechanical structures are presented, using an extended Picard-type iteration. Explicit discrete-mechanics approximations are taken as starting guess for the iteration. Iteration and necessary symbolic operations need to be performed only before time-stepping procedure starts. In a previous investigation, we demonstrated computational advantages for free vibrations of a hanging pendulum. In the present paper, we ﬁrst study forced non-linear vibrations of a tower-like mechanical structure, modeled by a standing pendulum with a non-linear restoring moment, due to harmonic excitation in primary parametric vertical resonance, and due to excitation recordings from a real earthquake. Our technique is realized in the symbolic computer languages Mathematica and Maple, and outcomes are successfully compared against the numerical time-integration tool NDSolve of Mathematica. For out method, substantially smaller computation times, smaller also than the real observation time, are found on a standard computer. We ﬁnally present the application to free vibrations of a hanging double pendulum. Excellent accuracy with respect to the exact solution is found for comparatively large observation periods.


Introduction
The need for novel computationally efficient time-integration schemes is increasingly invoked by online and real-time applications of non-linear dynamics of mechanical structures. Developments are rapidly proceeding to extended applications, e.g., in earthquake excited structures, concerning control of structural vibrations, see, e.g., [1,2], hybrid testing, see, e.g., [3], or non-linear real-time hybrid simulation, see, e.g., [4], for a force correction method. For the development of dynaimc structural models in linear non-linear earthquake engineering, as well as for suitable time-integration methods and non-linear control formulations, the reader is moreover referred exemplarily to [5][6][7][8][9][10]. In the following, a novel efficient time-integration scheme recently derived by the present authors, see [11], is applied for simulating non-linear forced and free vibrations of mechanical structures, and its advantages against standard numerical time-integration methods, such as NDSolve, available in Wolfram Mathematica [12], are exemplarily demonstrated. Our technique consists in an application of an extended Picard iteration to the time-integrated (global balance) form of the equations of motion, see [13] for the original Picard iteration and [14] for the structural equations of motion in integral form, i.e., for the global relations of balance of momentum. Concerning the applicability of this technique, we believe that it will be suitable as long as the derivation of the equations of motion has lead to a system of linear or non-linear ordinary differential equations of second order for the required degrees-of-freedom (DOFs), accompanied by the necessary number of initial conditions. In structural mechanics, such an initial-value system can be obtained by starting from the relations of (linear and/or angular) momentum for suitably modelled substructures (rigid bodies, deformable finite elements after discretization), taking into account constitutive relations, and using the usual reduction techniques to obtain a system with minimum number of DOFs. The system afterwards must be formally time-integrated to a system of first order. The time-integrated system then represents a system of so-called the global relations of balance of momentum, or first integrals. Often, when possible, it is more advisable to use the latter from the scratch, due to their wider applicability, since only the existence of integrals are to be required from a mathematical point of view, and no differentiability, or bounded integrands. Only mild mathematical assumptions hence must be made for allowing to substitute approximations into the integrands of these global balance forms of the equations of motion to obtain an improved solution, as this was originally suggested by Picard [13]. In the following, in extension of Picard [13], explicit discrete mechanics-type approximations, see Greenspan [15], or Runge-Kutta approximation methods, see [16], are used as starting guess in the Picard iteration, together with some advances that are nowadays offered by symbolic computation tools, e.g., in [12]. A main feature is that, thanks to the power of symbolic methods, our technique can be formulated so that the Picard iteration and the necessary symbolic operations need to be performed only once, before the time-stepping procedure starts. In our previous investigation [11], computational advantages have been demonstrated for large free vibrations of a hanging rigid pendulum, for which exact solutions do exist, and which is used in the literature as a benchmark example for comparison with numerical time integration methods. It was shown in [11] that a large number of free non-linear vibration cycles can be accurately computed, where the simulation remains in the close vicinity of the exact free phase-plane trajectory of the pendulum, due to the algorithmic satisfaction of the relations of global balance. Excellent accuracy with respect to the exact solution and competing numerical schemes was observed, also for large observation periods. In the present paper, we first study the application to the non-linear forced motions of a tower-like mechanical structure in the form of a standing (inverted) rigid pendulum with a non-linear restoring moment under a ground excitation of the earthquake-type and gravity. Structural parameters are taken from a real tower discussed in [17]. For small vibrations of the earthquake-excited pendulum with an inelastic spring considering the P-Delta effect, see [18]. For some other advanced applications of pendulum vibrations in the present context, see, e.g., [19,20]. In the following, we consider two cases of excitation: a harmonic excitation, which in vertical direction is in primary resonance, and a real earthquake excitation from recorded acceleration data. At a smaller scale, the model of a standing pendulum is also of interest for damping ground-excited vibrations of mechanical structures, see [21]. In the present contribution, we utilize results of the fourth-order Runge-Kutta method [16] as guess in the first step in the extended Picard iteration. Very similar time recordings of the non-linear structural response of our method are observed for our method on a standard computer, when compared to the standard initial value problem (IVP) solver NDSolve [12], but with substantially smaller computation times. Since computation times are also much smaller than the real observation periods of earthquake with a long duration, the technique should be of interest for real-time applications, such as in automatic control, particularly also because of the underlying symbolic formulation. In the last part of the paper, as a first step towards applications to non-linear multi-body dynamic systems, we present the application to free vibrations of a hanging double pendulum, for which also exact solutions do exist. In this example, two explicit discrete-mechanics type solutions originally suggested by Greenspan [15] are used and compared as starting guesses in the first Picard iteration. Symbolic computations are performed in Maple [22]. In order to show that the method can also be used conveniently for linear dynamic problems, we restrict to free vibrations, for which closed-form solutions do exist, and concentrate on accuracy aspects here. Excellent accuracy with respect to the exact solution is demonstrated, also for a comparatively large observation period, which proves that the method can be applied with high advantage for linear vibrations also.

Rigid-Body Model of an Earthquake Excited Tower-Like Structure
A fundamental non-linear model of a tower-like structure is shown in Figure 1. We regard the tower-like structure as an inverted pendulum with a non-linear spring at its hinge support, which is excited in both vertical and horizontal directions. We denote the total mass of the structure by m and its radius of inertia by i. The base accelerations are a H and a V . The mass center is located at a distance L from the support. We take a non-linear cubic restoring moment for the spring: M = −k L ϕ + k N ϕ 3 . The horizontal and vertical support reactions are called F H and F V . The unknown degree-of-freedom (DOF) in this 1DOF model is the angle ϕ counted from the vertical direction. (1) Derivatives with respect to time t are indicated by superimposed dots. Balance of linear momentum may be written as We may use the latter relation in order to eliminate the support reactions from balance of moments of momentum Equation (1) and end up with the following non-linear secondorder ordinary differential Equation (ODE): This relation can be regarded as balance of moment of momentum about the moving support, see ( [23], Chapter 7).With the frequency Ω 0 of the freely hanging linear pendulum when π − ϕ 1, we obtain For moderately large rotations, this non-linear ODE may be simplified by truncating the Taylor expansions for trigonometric terms as follows: The corresponding cubic non-linear ODE reads This equation will be treated as benchmark problem subsequently. Its linearized version has the form where the squared linear frequency of the pendulum with a linear restoring spring is given by The primary instability regime for vertical parametric resonance, see, e.g., ([17], Chapter 23) and ( [23], Chapter 10), extends from the following frequency For the harmonic excitation at this parametric resonance frequency, we obtain the following non-linear second-order ODE, assuming the horizontal ground excitation frequency also to be in resonance, i.e., to be equal to the vertical parametric resonance frequency: Amplitudes of horizontal and vertical earthquake accelerations are denoted as a H0 and a V0 , respectively. It is convenient to convert Equation (10) into the normal form , We now present a discussion of the computational method published in [11] and how to apply it in the present context. Our method is a time-stepping procedure; in the following, we consider time intervals of equal length T. The normal form of the balance equation in Equation (11) is brought into its integral form by integrating over the finite interval (time step) τ ∈ [0, T], the local time in the time interval being defined as τ = t − (n − 1)T, where n = 1, 2, . . . is the number of the time step under consideration, and t denotes the physical time, starting at the beginning of the observation period. We thus obtain the following formal representation of the angular velocity for any local time τ in the interval The reason for assigning f 2 formally with a superimposed tilde is explained later, after Equation (15). The angle of the tower from the vertical direction then follows as The generally inhomogeneous initial conditions at the beginning of the time interval under consideration are: The two integral Equations (12) and (13) are solved in an iterative manner, motivated by the work of Picard [13]. We however perform the first iterative step by substituting into the integral in Equation (12) the result obtained for ϕ from suitable explicit discretemechanics-type scheme. As a reference one, we here take the explicit fourth-order Runge-Kutta time-integration scheme, see [16]: An analogous result is obtained for ω(τ), which is not repeated for the sake of brevity. Having substituted Equation (15) into Equation (12), the integration in Equation (12) is not performed directly, but is carried out via a truncated Taylor series representation of the Runge-Kutta approximation for ϕ(τ) of Equation (15). The number of Taylor series terms may change from iteration step to iteration step. Thus, the tilde forf 2 in Equation (12) refers to truncated Taylor series of some prescribed order. The integral in Equation (13) is solved by substituting ω of Equation (12), where τ is replaced byτ. The necessary analytic operations, such as Taylor-series representations and integration, can be easily performed by means of symbolic computation. Since we approximate the relations of balance in their integrated (global balance) form, and not directly in their differential one and since we start the iteration using the well-established explicit fourth-order Runge-Kutta scheme, convergence is generally fast; for large free non-linear vibrations of the pendulum, see the detailed comparative study in [11]. The necessary analytic operations of the procedure, such as Taylor-series representations and integration, can be easily performed by means of symbolic computation. At the end of iteration, we obtain an analytic formula for the time evolution of ϕ(τ) and ω(τ) during the time interval τ ∈ [0, T]. Setting τ = T yields the values at the end of the time interval, ϕ(T) and ω(T). In a time stepping procedure, these values serve as initial values for the next time step. The analytic results obtained in the first time step after the end of iteration can be directly utilized, where τ now means local time in the next time step, starting again at the beginning of the latter.
A main advantage of this technique is that the Picard-type iteration and the analytic operations necessary for performing it must be done only once, for the first time interval, i.e., before the time-stepping procedure starts. Analytic expressions for ϕ and ω then have been obtained by means of symbolic computation, which can be stored and used in every time step. The considered numbers of Taylor terms, the number of iterations, and the length T of the time steps serve as parameters of the numerical problem, which allows adjusting the cost/accuracy ratio of the proposed computation for a specific problem at hand.

Oscillations under Harmonic Earthquake Excitation in Both Directions
For a realistic computation, we take structural parameters that we estimated from a special tower-like structure discussed in Chapter 23 of the monograph by Petersen and Werkle [17]: i = 8.66 m is the radius of inertia, L = 15 m is the distance of mass, mg = 3600 kN is the total gravity force, k L = 3 × 10 10 Nm is the linear stiffness coefficient, and we take k N = 0.9 × 10 10 Nm a non-linear stiffness coefficient. Moreover, we use a H0 = 9.81 m/s 2 as the horizontal acceleration amplitude, a V0 = 2.94 m/s 2 as the vertical acceleration amplitude, and T = 0.1 s is the considered time step length.
We perform the computations using the method described above (Picard-type iterations in Mathematica) and compare with the numerical standard function NDSolve of Mathematica, see [12]. The algorithm is described in Appendix A.1. The comparison is given in Figure 2. The results of both methods agree well. The proposed scheme has an advantage in the computation time: 0.20 s versus 0.94 s. In the derivation of the explicit formula used to plot Figure 2, in our method we implement three subsequent Picard-type iterations with 4, 6, and 9 Taylor terms, respectively, and compile the derived formula using the procedure Compile of Mathematica. Increasing the number of Picard iterations and the terms in the Taylor expansion would not change the relative error further in a visible manner. We refer the reader interested in technical details how the results are affected by the number of terms of Taylor expansion and by the time step length to Appendix A.2.

Numerical Results for Real Earthquakes
To perform the numerical computations, we again use our method and the standard initial value problem solver NDSolve of Mathematica, see [12]. We exemplarily use a ground acceleration record from the Kobe earthquake, considered to be one of the most devastating and costly natural disasters in recent history [24]. The Kobe earthquake data, originally given as equally spaced acceleration values, on which our computations are based, are shown in Figure 3. For the NDSolve computations, earthquake data are interpolated by means of the standard procedure Interpolation of Mathematica (the interpolation order is set to 0, which means a piecewise-constant function at the output), see also Figure 3. The values from the earthquake data are used as constants in the corresponding time steps. For simplicity sake, the vertical excitation is taken proportional to the horizontal one, a V = 0.3a H . The set of Kobe earthquake data has the time step T = 0.01 s. The number of time steps used in the computations is 4090. We use two iterations with 7 Taylor terms each and with standard machine precision in the computations and compile the resulting symbolic formulas using the procedure Compile of Mathematica. The results of the numerical integration in the case of real earthquake record in Figure 3 are plotted in Figure 4, which shows the comparison between the results of our method and the NDSolve procedure of Mathematica. We observe well agreement between the results, where the computation time is substantially shorter for the proposed method. An analogous behavior was found when using records for other disastrous earthquakes, such as the El Centro earthquake and the Bam earthquake. We note, however, that our choice of interpolating earthquake data as piecewise-constant, which was motivated by the practical requirement not to add further functional assumptions to the point-wise measured data, may impose significant performance penalty to NDSolve. Namely, it forces this algorithm to use extremely short steps at the jumps. In contrast, for the present Picardtype method the latter jumps do not result in problems, because, as mentioned already above, differentiability of the integrands is not required, only existence of the integrals.

Double Pendulum
In this section, we investigate free vibrations of the double pendulum using the above presented method and compare the results with the analytic solutions. The scheme of the double pendulum is given in Figure 5.
The system of equations reads where the dimensionless time is introduced as The problem is linear and admits the exact solution in the form in the special case with the following initial conditions For the double pendulum we first apply a discrete-mechanics-type scheme by Greenspan [11,15], subsequently denoted as Greenspan I, as the first guess in the Picardtype iteration. The calculation technique consists of the following steps: The last two formulas are then integrated iteratively four times with the integral form of the balance of momentum and with the integral kinematic relation where In this section, we use the computer algebra system Maple [22] for symbolic computations and plotting. In Figure 6, we draw the relative single-step error with the Greenspan I scheme for the double pendulum. The terms of the approximate solution with k iterations are the same up to the 2k + 2 terms of the Taylor expansion of the exact solution. Even for large time steps is the error rather small. We also apply the second discrete-mechanics-type Greenspan scheme (Greenspan II) [11,15] given by , with the mid-point formula for the angle The result is then iteratively integrated four times with the balance relation and with the kinematic relation.
In Figure 7 we present the relative errors in the angles for the Greenspan II scheme. We see in the plots that the relative errors are again small even for larger time steps. The terms of the approximate solution with k iterations are the same up to the 2k + 2 terms of the Taylor expansion of the exact solution. There are further terms in the approximate solution, which doe not exactly correspond to the exact solution. The errors of four iterations are presented; the results appear to be converging to the exact solution as expected. Our last task is to present the relative errors in a multistep procedure composed from the Greenspan II scheme with four iterations. The results are given in Figure 8 for a large number of oscillations corresponding to n =10,000 time steps.

Conclusions
In the above exemplary study, we showed the application of our discrete-mechanicstype time-integration technique with extended Picard-type iterations, first published in [11], to the example of an earthquake-excited tower-like structure, modeled as a rigid non-linear pendulum. The corresponding vibration problems were also solved using the standard ND-Solve method Wolfram Mathematica [12], and the outcomes were compared successfully to our method, where our proposed approach required considerable less computation time. Free vibrations of the double pendulum were studied afterwards in order to demonstrate that the method can be applied with advantage under linear conditions also, leading to a high accuracy even for comparatively large time-steps. Since our explicit time-integration method is based on analytic expressions, which can be obtained by symbolic computation in a straight-forward manner, it is hoped that it will find interest in various scientific engineering applications, e.g., in earthquake engineering, in vibration control, and for hybrid testing.