A locally defined time-marching technique for structural dynamics

Abstract. In this work, a new time marching procedure is proposed for structural dynamics analyses. In this novel technique, time integration parameters are locally defined and different values may be attributed to each structural element of the model. In addition, the time integrators are evaluated according to the properties of the elements, and the user may select in which structural elements numerical dissipation will be introduced. Since the integration parameters are locally defined as function of the structural element itself, the time marching technique adapts according to the model, providing enhanced accuracy. The method is very simple to implement and it stands as an efficient, direct, single-step procedure. It is second order accurate, unconditionally stable, truly self-starting and it allows highly controllable algorithm dissipation in the higher modes. Numerical results are presented along the paper, illustrating the good performance of the new technique.


Introduction
The development of direct time integration schemes for both linear and nonlinear structural dynamics analysis has been the subject of several research works in the last decades.Since the response of interest in structural dynamics applications is usually located in low frequencies modes [1], a desirable and important property of time marching procedures is the introduction of numerical damping which suppresses the participation of artificial high frequencies (derived from spacial discretization) without affecting the low frequencies modes.In this sense, many dissipative and non-dissipative algorithms are reported in the literature [1][2][3][4][5][6][7][8][9][10].The classical Newmark method [3] is widely used in structural dynamics problems, however it is first order accurate when numerical damping is introduced.The Bathe method [7] is dissipative and remains second order accurate, nevertheless it is a two step method (i.e., it requires two systems of equations to be dealt with within each time step), and presents large amplitude decay in long time analysis.Thus, following the work presented by Soares [11,12], a new methodology, in which the user may control the introduction of numerical dissipation at the element level, is proposed.The basic equations of the new time marching procedure are discussed here, and the obtained numerical results of the proposed approach are compared to those provided by the Newmark and the Bathe methods.

Governing equations and time integration scheme
The governing system of equations for structural dynamics problems may be written as [13]: where M, C and K are mass, damping and stiffness matrices, respectively; F(t) is the applied load vector and U(t), U(t) and Ü(t) are the displacement, velocity and acceleration vectors, respectively.The initial conditions are given by: U(0) = U 0 and U(0) = U0 ; where U 0 stands for the initial displacement vector and U0 stands for the initial velocity vector.Considering a constant time step ∆t and the integration of Equation 1 with respect to time, at the element level (represented with subscript e), the following expression is obtained: The integrals in the left-hand side of Equation 2 may be approximated as [11]: where superscripts n and n + 1 indicate the time step of the variables.Equation 3c is the chosen approximation of the integral of U e (t) [11].Both α n e and γ n e are the parameters of new methodology, calculated for each element at each time step.The displacement U n+1 is defined as: Considering approximations in 3 and 4, Equation 2 may be rewritten as: (5) Thus, after properly assemblage, velocities may be computed following Equation 5 and displacements may be evaluated following Equation 4, establishing the proposed time marching procedure.Once the new methodology is based only on velocity and displacement relations, no computation of accelerations is required and therefore it is a self-starting method.Integration parameters γ n e and α n e are locally calculated and different values may be attributed to each element, at each time step of the analysis.Since only linear analysis will be considered here, these parameters remain the same along time (γ n e = γ e and α n e = α e ).Such parameters are computed based on the properties of the elements and the user may select in which element numerical dissipation will be introduced.In order to achieve this, each structural element will be provided with a numerical dissipation property, called a e .Whenever numerical dissipation is necessary, it must be adopted a e > 0, otherwise a e = 0.Then, if a e = 0, γ e and α e are computed as follows [11,12]: where ω max e stands for an approximation for the maximum natural frequency of the element.Otherwise, if a e > 0, γ e and α e are given by [11,12]: tanh a e ω max e ∆t (7a) introducing numerical dissipation into the analysis, at the element level.Thus, the intensity of the numerical dissipation is locally controlled by the user, according to the given a e value (the greater the a e , the greater the numerical damping feature).
The expressions described by Equations ( 6) and ( 7) are deduced in order to keep the proposed technique unconditionally stable, second order accurate, and with optimal dissipative features (following Equation (7b), minimal values of spectral radii occur at maximal sampling frequencies).In addition, when Equations ( 6) are followed, describing a non dissipative methodology, the new technique can be proven to be always more accurate than the standard Trapezoidal Rule (in fact, the Trapezoidal Rule is reproduced by the new approach when ω max e tends to infinity, following Equation (6a)).

Numerical application
This numerical application has been extracted from [14] and consists of a simple three degrees of freedom spring model, as shown in Figure 1.Since node 1 is subjected to a prescribed displacement (u 1 = sin ω p t), it is possible to write the effective governing equations of the model as: where ω p = 1.2, k 1 = 10 7 , k 2 = 1, m 2 = 1, m 3 = 1, and null initial conditions are considered.The adopted time step was ∆t = 0.2816s, and a long period of analysis was considered (i.e., 100s), to better illustrate the period elongation and amplitude decay errors of the computed results.Taking into account the new approach, a e = 1 was adopted for the first element of the model, and a e = 0 was considered for the second.The integration parameters γ e and α e are then automatically computed according to Equations ( 6) and (7), and the obtained values are described in Table 1.https://doi.org/10.1051/matecconf/201821117004VETOMAC XIV Table 1.Numerical dissipation property adopted for each element and respective integration parameters.
Element a e γ e α e 1 st 1.0 2 1 2 nd 0.0 0.0461 0.9539 Despite its simplicity, this model may be used to reproduce a much more complex structural system.The left high stiffness spring may represent an almost rigid connection while the left flexible spring represents the flexible parts of a complex system [14].Figures 2 and  3 give the calculated solutions for velocities at node 2 and 3, respectively.Since the Bathe method is a two step method (i.e., it requires two systems of equations to be dealt with within each time step), the adopted time step for Bathe method was twice larger than the others (∆t = 0.5632s), in order to evenly compare the computed results in terms of computational costs.In addition, a reference response, obtained by a mode superposition solution using the lowest frequency and static correction [14,15], is presented.
Relative error results are evaluated as follows: where u(t j ) stands for the computed solution, u A (t j ) stands for the reference answer and N t corresponds to the total number of time steps in the analysis.Error results are presented in Table 2.As one may observe, the proposed technique presented better accuracy, fast dissipation of spurious high frequencies (Figure 2(b)) and lower period elongation and amplitude decay errors (Figure 3(c)).It is important to observe that this enhanced performance is obtained considering the same computational effort of the adopted comparing methodologies, highlighting the great effectiveness of the proposed approach.

Conclusions
In this work, a novel time marching procedure for structural dynamics, following the work presented by Soares [11,12] is discussed.In section 3, numerical results are presented, illustrating the better performance of the new method.The main advantages of the proposed technique are: (i) it is simple; (ii) it demands low computational costs, since it stands as a single step procedure, considering systems of equations with regular (i.e., not expanded) dimensions; (iii) the user may control numerical dissipation at element level; (iv) it is second order accurate; (v) it is unconditionally stable; (vi) it is truly self-starting; (vii) integration parameters are locally evaluated, based on the properties of the element, enabling a more effective procedure.

Table 2 .
Relative errors for velocities at node 2 and 3.