The Direct Integration Method with Virtual Initial Conditions on the Free and Forced Vibration of a System with Hysteretic Damping

: The energy dissipated in hysteretic damping is independent of cyclic frequency, which agrees with the experimental results of energy dispersion of many materials under cyclic loading. Despite these desirable properties, hysteretic damping su ﬀ ers from a notable drawback; the direct integration solution of the equation of motion is divergent. In this paper, a virtual initial condition is proposed to address the instability of the solution. The new method develops virtual initial conditions associated to the real initial conditions, which make the direct integration solution eliminate the divergence term of the complementary solution and converge to the exact solution. The stability and accuracy of the proposed direct integration method are demonstrated to be e ﬀ ective in solving the divergence problem by the comparison of the numerical and theoretical solutions on the free and forced vibration of a system with hysteretic damping.


Introduction
Damping is the natural property of a vibration system, which dissipates vibration energy. The widely used models to represent damping in a system are the viscous model and the hysteretic model [1][2][3]. Viscous damping, whose force is proportional to the velocity, is the most popular model due to the mathematical convenience. However, it implies a dissipated energy per a vibration cycle depending on the frequency. Hysteretic damping, whose force is proportional to the displacement and in phase with the velocity, leads to a dissipated energy per cycle which is independent of the frequency [4,5]. Experiments on structural materials indicate that the energy dissipation is independent of the cyclic frequency [6][7][8]. In those cases, the hysteretic damping is more precise. However, there are some drawbacks in the solution of the equation of motion with hysteretic damping for both the frequency domain analysis and the direct integration method in the time domain, which hinders the wider application of hysteretic damping models.
The frequency domain approach solves the hysteretic damping system by using the Fourier transform for the equation of motion [9][10][11], but it has two main drawbacks [12,13]. (1) Part of transient response may be neglected which may underestimate the response at the initial time and (2) in theory, it is only applicable for the linear elastic system due to the application of the superposition principle, so it cannot be used to analyze the true nonlinear response. For the approximate solution of a nonlinear system, Idriss and Seed [14] proposed an equivalent linearization method, which was widely used in the nonlinear seismic response of soil site and soil-structure interaction [15][16][17][18].
In a time-domain analysis, the numerical method is widely used because of its high efficiency, and because it can analyze the nonlinear problems as well [19][20][21][22][23]. However, the result from the Appl. Sci. 2019, 9, 3707 2 of 12 numerical method is divergent when the direct integration method is used to solve the hysteretic damping systems, even if the integration procedure is unconditionally stable [24]. To develop the hysteretic damping in the time domain, a common method is to transfer the hysteretic damping into the approximate viscous damping. Henwood [25] used the forcing frequency to replace the natural frequency. Thus, the hysteretic damping matrix can be transferred to a viscous damping matrix, which is independent of frequency. Chen [26] proposed an integral-differential equation (IDE) in the time domain for the free vibration of a single degree-of-freedom (SDOF) system by using viscous damping depending on the frequency to equivalent hysteretic damping. Zhou [27] combined the standard linear constitutive models to approximate the hysteretic damping constitutive model in the time domain, and then applied the time-step integral calculation to hysteretic damping model.
However, until now, the improvement in the direct integration of a hysteretic damping system is still very slow, and the divergence of the results of the direct integral method has not been solved. In theory, there are two kinds of unstable solutions for the direct integration method. One is caused by the instability of the numerical method, and the other is caused by the divergent solution of the equation itself. The instability of the direct integration method of the hysteretic damping model is caused by the latter one since the pair of eigenvalues of the SDOF hysteretic damping are complex numbers and opposite to each other. The complementary solution of the equation is a linear combination of the two eigenvalue exponential functions. Therefore, the eigenvalue with a positive real part will cause a divergent solution, which is contrary to the natural and physical phenomenon. In the theoretical solution, the divergent term can be arbitrarily omitted, but it is kept in the direct integration procedure. Recently, Sun [28] proposed a time-domain numerical method based on the interpolation of excitation, which is a semi-analytic direct integration method, by using the theoretical solution in each step.
Therefore, the key to the stability of the direct integration method of a hysteretic damping system is to eliminate the divergent term in the recursive process. In this paper, we propose to use virtual initial conditions to address the instability of the solution. The proposed method constructs virtual initial conditions by the real part of the displacement and velocity of the previous step, as well as the load condition, so that the direct integration solution converges to the exact solution. Except for the initial conditions, all the procedures of the proposed method are the same as the direct integration method for the governing equation of motion with viscous damping. Then, the stability and accuracy of the proposed direct integration method are demonstrated by the comparison of numerical and theoretical solutions on the free and forced vibration of the SDOF system with hysteretic damping.

Equation of Motion and the Direct Integration Method
The governing equation of motion of a single-degree-of-freedom (SDOF) system with hysteretic damping in forced vibration can be expressed as [29] m ..
where m and k are the mass and stiffness, respectively, η is the hysteretic damping factor, i = ..
u and u are the acceleration and displacement of the SDOF, respectively, and f (t) is the applied force. This equation is to be solved for the displacement, u(t), subjected to the initial conditions as where u 0 and v 0 are the displacement and velocity at the time t = 0. For the direct integration method, let the time axis be divided into small intervals of equal length ∆t, and let t n = n∆t(n = 0, 1, 2, . . . . N) be the division times. The response of displacement, velocity, and acceleration at the discrete time instants t n are u n , . u n , and .. u n . The applied force f (t) at t n is f n .
So far, several direct integration methods have been proposed to obtain the solution at time t n+1 assuming that the solutions at times t i (i = 0, 1, 2, . . . . n) are known. The Newmark method is used to Appl. Sci. 2019, 9, 3707 3 of 12 investigate the stability of the solution of the hysteretic damped equation in the following illustration. In the Newmark method, the velocity and acceleration at time t n+1 can be expressed as ..
in which parameter γ and β are parameters that can be determined to obtain integration accuracy and stability. Substituting Equation (4) into Equation (1) gives: ..
u n ]m. Solve the displacement, u n+1 , by Equation (5), then . u n+1 may be determined from Equation (3) and the acceleration may be determined from the following expression deduced from Equation (1): It is well known that the Newmark method is unconditionally stable for γ = 1 2 and β = 1 4 , which is called the average acceleration method. However, the results obtained by the average acceleration method will blow up for the hysteretic damping system, giving meaningless results. The reason is that the solution will converge to the divergent solution which is caused by the divergent member of the complementary solution of Equation (1) if the real-valued initial condition is used directly in the Newmark method. Therefore, in this paper, a modified initial condition in each time step is proposed to address the unstable phenomenon for direct integration.

The Virtual Initial Condition for Free Vibration
For a system in free vibration, the right-hand member of Equation (1) vanishes and gives m ..
The solution for the homogeneous equation is in which C is a constant and λ is eigenvalue. Substituting Equation (8) into Equation (7) leads to the algebraic equation: The solutions of Equation (9) are . The two eigenvalues are opposite numbers to each other. The general solution of Equation (7) is Appl. Sci. 2019, 9, 3707 4 of 12 The constant C 1 and C 2 in terms of the real-value initial conditions of Equation (2) may be determined from Since the real part of s is a positive quantity, there would be a growing response along the time for the right-hand first member C 1 e st , which is a divergent term which is contrary to the physical rule and should be omitted. Therefore, by introducing C 1 = 0 the correct general solution can be expressed as where C is a complex constant. Applying the initial condition of Equation (2) gives In the analytical solution, the divergent member can be omitted to obtain the exact solution of Equation (14). However, the direct integration solution converges to Equation (12). Therefore, it is crucial to make Equation (12) equal to Equation (14). In this paper, the virtual initial conditions u v1 0 and v v1 0 are proposed to achieve this goal. The total initial conditions are the sum of the real initial conditions and corresponding virtual initial conditions, that is Substituting Equation (16) into Equation (13), and making Equation (12) equal to Equation (14), gives Solving the two algebraic Equations (17) and (18) leads to It can be seen that the virtual initial conditions are purely imaginary. The usual real initial conditions are It can be observed that the complete initial conditions are complex, but only the real parts can be observed and measured. The corresponding unobserved imaginary part should be the indivisible member of the initial condition. After introducing the virtual initial conditions, the general solution of Equation (7) is naturally exact, the real part of which is the observable solution and free decay.

The Virtual Initial Conditions for Harmonic Force
If f (t) is a harmonic force with forcing frequency θ, written as Ae iθt , the particular solution of Equation (1) is Appl. Sci. 2019, 9, 3707 The complementary solution is the free vibration response given by Equation (12). The complete solution is the sum of the complementary and particular solutions.
The constants D and E are determined by the initial conditions. To illustrate conveniently, the virtual initial conditions associated with the real zero initial condition are deduced first. On applying the zero initial condition with u(0) = 0 and . u(0) = 0 to Equation (23), the constants D and E are deduced.
Similar to the free vibration, the right-hand first term De st is the divergent term and contrary to the physical rule, which should be omitted. The complete solution that is satisfied with physical meaning is On making use of u(0) = 0 and . u(0) = 0 to Equation (25), the constant C f is deduced as where Re(X) is the real part of X, and Im(X) is the imaginary part of X.
The direct integration solution converges to Equation (23), therefore, let Equation (23) be equal to Equation (25) by making use of the virtual initial condition u v2 0 and v v2 0 . Subjected to initial displacement u v2 0 and initial velocity v v2 0 , the constants D and E in Equation (23) are determined from Solving the two algebraic Equations (27) and (28) leads to On making use of the identity between Equations (23) and (25), one obtains Then, it leads to The virtual initial conditions u v2 0 and v v2 0 are also purely imaginary. The member C f e −st in Equation (25) is free adjoint vibration which is independent of the real initial conditions but depends on the load. u v2 0 and v v2 0 are artificial load-dependent initial conditions which are developed to make the direct integration solution converges to the exact solution.
Equations (31) and (32) are deduced by the assumption that the system is initially at rest. For a system with non-zero initial conditions, the virtual initial conditions defined by Equations (31) Appl. Sci. 2019, 9, 3707 6 of 12 and (32) should be augmented by the addition of the free vibrational solution defined by Equation (19) and (20). Then the total initial conditions of harmonic responses are

The Virtual Initial Conditions for Arbitrary Force
Assuming the SDOF system with hysteretic damping is subjected to real-valued force g 1 (t), which varies arbitrarily with time, it is necessary to establish a dual force in the time domain analysis [30]. Then, the governing equation of motion can be written as m ..
in which g 1 (t) is the observable force, and g 2 (t) is the corresponding dual force. By using discrete Fourier transforms (DFT), g 1 (t) can be expanded to the sum of a series harmonic function, that is where θ j = 2π j N∆t (j = 0, 1, 2, . . . , N/2−1,N/2), ∆t is time intervals of discrete, Fourier coefficients A 0 , A 1 , A 2 , . . . , A N/2 and B 1 , B 2 , . . . , B N/2−1 can be determined from where t i = l∆t (l = 0, 1, 2, . . . , N−1) is the discrete value of t, g 1 (t l ) is the discrete value of g 1 (t). Then dual loading g 2 (t) can be expressed as So, Equation (39) denotes that arbitrary force can be expressed as the sum of a series of harmonic forces, so the complete solution is the sum of a series of harmonic forces.
The right-hand first term D j e st is a divergent term and contrary to the physical rule and should be omitted, so the complete solution which is satisfied with physical meaning is  1, 2, . . . , N/2−1), and X N/2 = A N/2 2(−θ 2 N/2 m+(1+iη)k) .
Similar to the harmonic force, the virtual initial condition associated with each term can be determined from Equations (31) and (32) for the zero initial condition. Therefore, the virtual initial condition for arbitrary force can be expressed as Then the total initial conditions for a system with non-zero initial conditions are Algorithm 1 summarizes the above-described procedure, which can be implanted into the computer. The procedure can also be used to calculate the free vibration and harmonic force with slight modification. Let f n = u v3 n = v v3 n = 0 for free vibration. The harmonic force is the particular example of the arbitrary force, which only keeps one term with the with forcing frequency θto calculate u v3 n and v v3 n . It is noted that step 2.4 is used to calculate the virtual initial condition for each time step n to eliminate the effect of the round-off and improve the accuracy. Then the stability is dependent on the β, γ. Algorithm 1 Newmark method with the virtual initial condition for arbitrary force.
Step 1.0 Initial calculation 1.1 DFT for g 1 (t) to calculate coefficients A j and B j by Equations (36) and (37).
Step 3.0 Repetition for the next time step. Replace n by n + 1 and implement step 2.1 to 2.6 for next time step.

Free Vibration
This example case study analyzes an SDOF with free vibration, which has m = 1, k = 4π 2 , and η = 0.1, considering the free vibration with the initial condition: Appl. Sci. 2019, 9, 3707 8 of 12 From Equation (14), the correct theoretical solution is However, the wrong solution by Equation (12) with only the real initial conditions is In this example, the problem was solved by the Newmark method with γ = 1 2 , β = 1 4 . The numerical results obtained using ∆t = 0.01s are compared with the theoretical solution, as shown in Figure 1, in which NMV is numerical result obtained by algorithm 1, NMR is numerical result obtained by the Newmark method with real-valued initial conditions. From Equation (14), the correct theoretical solution is Equation (16) converges to the correct theoretical solution of Equation (47). Therefore, the Newmark method is stable, but the solution is unstable for the real initial conditions. After applying the proposed virtual initial condition, the divergent term is eliminated to make the solution stable.

Harmonic Vibration
This example case study analyzed an SDOF subjected to harmonic force, i t e θ , with / The solution of the steady-state response is 1.571 (0.027 0.0028 ) The problem was also solved by the Newmark method with  It can be seen from Figure 1 that if directly using the real initial conditions u 0 = 1 and v 0 = 0 by the Newmark method, the displacement blows up to give meaningless results which converge to the wrong solution of Equation (48). However, the displacement with initial conditions u t 0 and v t 0 by Equation (16) converges to the correct theoretical solution of Equation (47). Therefore, the Newmark method is stable, but the solution is unstable for the real initial conditions. After applying the proposed virtual initial condition, the divergent term is eliminated to make the solution stable.

Harmonic Vibration
This example case study analyzed an SDOF subjected to harmonic force, e iθt , with θ/ω = 0.25, The solution of the steady-state response is u = (0.027 − 0.0028i)e 1.571it .
The problem was also solved by the Newmark method with γ = 1 2 , β = 1 4 . The numerical results obtained using ∆t = 0.01s are compared with the theoretical solution, as shown in Figure 2, as well as the solution of steady-state response. Obviously, the displacement blows up again for initial conditions u 0 = 0 and v 0 = 0, and it converges to the wrong solution of Equation (50). The result is stable and converges to the correct theoretical solution of Equation (49) for initial conditions u t 0 and v t 0 . Therefore, for forced vibration, the unstable numerical solution is induced by the free adjoint vibration, which includes the divergent term. The virtual initial condition removed the divergent term so that the numerical solution is the same as the theoretical solution. v . Therefore, for forced vibration, the unstable numerical solution is induced by the free adjoint vibration, which includes the divergent term. The virtual initial condition removed the divergent term so that the numerical solution is the same as the theoretical solution.
The numerical results of NMV are the total response which includes transient vibration and steady-state vibration. The transient response will decay exponentially with time, and only the steady-state vibration will remain. However, the peak displacement of the correct theoretical solution by Equation (49) is about 1.6 times than that of steady-state response at the initial stage. Since the solution in the frequency domain is only the steady-state vibration, it represents that neglecting the effects of transient response would lead to unreliable results in engineering design.

Seismic Excitation
This example case study analyzed an SDOF with m = 1, k = 4π 2 , and η = 0.1 subjected to seismic excitation with ground acceleration ( ) . The equivalent real-valued force g1(t) is In this example, acceleration ( ) g u t  was the ground motion recorded at a site in Whittier Narrows, California during the Whittier earthquake of 1 October 1987. The force g1(t) is shown in Figure 3.   The numerical results of NMV are the total response which includes transient vibration and steady-state vibration. The transient response will decay exponentially with time, and only the steady-state vibration will remain. However, the peak displacement of the correct theoretical solution by Equation (49) is about 1.6 times than that of steady-state response at the initial stage. Since the solution in the frequency domain is only the steady-state vibration, it represents that neglecting the effects of transient response would lead to unreliable results in engineering design.

Seismic Excitation
This example case study analyzed an SDOF with m = 1, k = 4π 2 , and η = 0.1 subjected to seismic excitation with ground acceleration .. u g (t). The equivalent real-valued force g 1 (t) is In this example, acceleration .. u g (t) was the ground motion recorded at a site in Whittier Narrows, California during the Whittier earthquake of 1 October 1987. The force g 1 (t) is shown in Figure 3. v . Therefore, for forced vibration, the unstable numerical solution is induced by the free adjoint vibration, which includes the divergent term. The virtual initial condition removed the divergent term so that the numerical solution is the same as the theoretical solution.
The numerical results of NMV are the total response which includes transient vibration and steady-state vibration. The transient response will decay exponentially with time, and only the steady-state vibration will remain. However, the peak displacement of the correct theoretical solution by Equation (49) is about 1.6 times than that of steady-state response at the initial stage. Since the solution in the frequency domain is only the steady-state vibration, it represents that neglecting the effects of transient response would lead to unreliable results in engineering design.

Seismic Excitation
This example case study analyzed an SDOF with m = 1, k = 4π 2 , and η = 0.1 subjected to seismic excitation with ground acceleration ( ) . The equivalent real-valued force g1(t) is In this example, acceleration ( ) g u t  was the ground motion recorded at a site in Whittier Narrows, California during the Whittier earthquake of 1 October 1987. The force g1(t) is shown in Figure 3.     Figure 4b shows the Newmark method with the real-valued initial conditions (NMR) and the wrong solution by Equation (40) with the real initial conditions. The Newmark method's parameters are γ = 1 2 , β = 1 4 , and the integration time interval is ∆t = 0.02s.
the Newmark method with the real-valued initial conditions (NMR) and the wrong solution by Equation (40) with the real initial conditions. The Newmark method's parameters are  Figure 4a,b further confirm that the divergence of the numerical solution is caused by the divergent term of complementary solution. The Newmark method would converge to the wrong theoretical solution for initial conditions 0 0 u = and 0 0 v = . However, after introducing the virtual initial conditions, the Newmark method would converge to the correct theoretical solution. The difference between the NMV and frequency domain method is in the initial stage. Therefore, the frequency domain method would obtain reasonable results if the response is relatively small at the beginning of vibration.

Conclusions
In this paper, a direct integration virtual initial conditions method is developed for the numerical solution of the free or forced vibration of a hysteretic damped system. Based on extensive theoretic derivations and numerical analysis, the following conclusions can be drawn: (1) For free or forced vibration of a hysteretic damped system, the stability of numerical methods is the same as that of a viscously damped system. The divergence of numerical results is caused by the divergent term of the complementary solution if using the real-valued initial conditions. The virtual initial conditions can remove the divergent term, and make the numerical solution converge to the exact theoretical solution. (2) The virtual initial conditions are purely imaginary. For free vibration, the virtual initial conditions depend on the real-valued initial conditions, and for forced vibration, the virtual initial conditions depend on the amplitude of the force. (3) The solution in the frequency domain can only obtain the steady-state vibration. However, the proposed direct integration method can accurately calculate the transient response, which results in a reasonable estimation for whole vibration history and is, thus, recommended for practical applications.
In the future, the direct integration virtual initial conditions method is recommended to be further developed for the multi-degree-of-freedom system with hysteretic damping and non-linear dynamic response.
Author Contributions: Conceptualization, D.P.; methodology, D.P. and X.F.; investigation, X.F. and W.Q.; writing-original draft preparation, X.F.; writing-review and editing, D.P.  Figure 4a,b further confirm that the divergence of the numerical solution is caused by the divergent term of complementary solution. The Newmark method would converge to the wrong theoretical solution for initial conditions u 0 = 0 and v 0 = 0. However, after introducing the virtual initial conditions, the Newmark method would converge to the correct theoretical solution. The difference between the NMV and frequency domain method is in the initial stage. Therefore, the frequency domain method would obtain reasonable results if the response is relatively small at the beginning of vibration.

Conclusions
In this paper, a direct integration virtual initial conditions method is developed for the numerical solution of the free or forced vibration of a hysteretic damped system. Based on extensive theoretic derivations and numerical analysis, the following conclusions can be drawn: (1) For free or forced vibration of a hysteretic damped system, the stability of numerical methods is the same as that of a viscously damped system. The divergence of numerical results is caused by the divergent term of the complementary solution if using the real-valued initial conditions. The virtual initial conditions can remove the divergent term, and make the numerical solution converge to the exact theoretical solution. (2) The virtual initial conditions are purely imaginary. For free vibration, the virtual initial conditions depend on the real-valued initial conditions, and for forced vibration, the virtual initial conditions depend on the amplitude of the force. (3) The solution in the frequency domain can only obtain the steady-state vibration. However, the proposed direct integration method can accurately calculate the transient response, which results in a reasonable estimation for whole vibration history and is, thus, recommended for practical applications.
In the future, the direct integration virtual initial conditions method is recommended to be further developed for the multi-degree-of-freedom system with hysteretic damping and non-linear dynamic response.