Abstract

Real-time substructure testing (RST) algorithm is a newly developed integration method for real-time hybrid simulation (RTHS) which has structure-dependent and explicit formulations for both displacement and velocity. The most favourable characteristics of the RST algorithm is unconditionally stable for linear and no iterations are needed. In order to fully evaluate the performance of the RST method in solving dynamic problems for nonlinear systems, stability, numerical dispersion, energy dissipation, and overshooting properties are discussed. Stability analysis shows that the RST method is only conditionally stable when applied to nonlinear systems. The upper stability limit increases for stiffness-softening systems with an increasing value of the instantaneous degree of nonlinearity while decreases for stiffness-hardening systems when the instantaneous degree of nonlinearity becomes larger. Meanwhile, the initial damping ratio of the system has a negative impact on the upper stability limit especially for instantaneous stiffness softening systems, and a larger value of the damping ratio will significantly decrease the upper stability limit of the RST method. It is shown in the accuracy analysis that the RST method has relatively smaller period errors and numerical damping ratios for nonlinear systems when compared with other two well-developed algorithms. Three simplified engineering cases are presented to investigate the dynamic performance of the RST method, and the numerical results indicate that this method has a more desirable accuracy than other methods in solving dynamic problems for both linear and nonliner systems.

1. Introduction

Real-time hybrid simulation (RTHS) is an advanced technique for assessing the seismic behaviours of structures, especially large-scale and ultralimit complex systems [14]. Similar to other hybrid simulation techniques, the original system in a RTHS is divided into two components, including numerical substructures and experimental substructures. The displacement-coordination and force-equilibrium are two principles for the coupling between two components at the interfaces [5]. Both the restoring force calculated from numerical substructures and the resisting force measured from experimental substructures at each time step are fed back to an integration algorithm. Then, the subsequent time-step target displacement applied to the experimental substructure is calculated by solving the step-by-step equation of motion under external excitations and the reacting forces measured from experimental substructures [6]. “Real-time” is a significant characteristic that makes RTHS different from other hybrid simulations, and the corresponding applications include “real-time calculating of numerical substructures”, “real-time loadings on experiment substructures,” and “real-time transferring the interactions between these substructures.” So far, a few works have focused on the latter two applications. However, with an increasing size and complexity of structural systems, “real-time calculating” becomes more important and the requirement of numerical integration algorithm is higher.

In RTHS, integration algorithms are effective methods to solve discrete equations of motion for structural dynamics [7]. Integration algorithms are usually classified as explicit and implicit based on the expressions of the displacement and the velocity. In a pseudodynamic testing that quasi-statically imposes the target displacement on the experimental substructures, an integration algorithm is explicit when the displacement for the next time step can be obtained from the response quantities (displacement, velocity, and acceleration) at the current and previous time steps, and otherwise it is implicit [8]. However, the definitions of explicit and implicit algorithms have changed in RTHS, as the loading rate becomes faster or real-time and the velocity of experimental substructures should be feedback for integration. Therefore, an integration algorithm is explicit when both the displacement and velocity for the next time step can be determined as a function of the response quantities at the current and previous time steps. Although various integration algorithms have been used to conduct RTHS, explicit integration algorithms [918] are generally preferred over implicit ones [1923] as they do not involve nonlinear iterations and show a higher calculation efficiency. However, conditional stability limits the application of conventional explicit algorithms in solving structural dynamic problems because it may lead to an extremely small step size for large-number multidegree-of-freedom (MDOF) systems which have high frequency modes.

More and more researchers focus on developing new types of integration algorithms, which can achieve both explicit formulation and unconditional stability at the same time. Chang [24] developed an explicit integration algorithm which is referred to as Chang explicit method (CEM) for pseudodynamic testing, where the expression of displacement is an explicit form. However, when used in RTHS, CEM becomes semiexplicit since the determination of the velocity for the next time step requires knowledge of the unknown acceleration at the next time step. The major difference between CEM and other algorithms lies in the coefficients of difference equations introduced in displacement increment which are functions of initial structural properties and time step rather than constants, making it a structure-dependent algorithm. Chen and Ricles [12] proposed another type of structure-dependent integration algorithm named CR method (CRM) based on discrete control theory. Different from CEM, the expressions for displacement and velocity are both in explicit forms and the coefficients for both displacement increment and velocity increment are functions of initial structural properties and time step. It is noted that CEM and CRM are unconditionally stable for linear and instantaneous stiffness softening systems and possess second order accuracy [25, 26]. Kolay and Ricles [14] introduced a family of explicit direct integration algorithms using discrete control theory, which is named as explicit KR-α algorithm. This method is unconditionally stable for linear and instantaneous stiffness softening systems, and it can control the amount of numerical damping by a single parameter. However, all these favourable attributes are achieved at the expense of degraded accuracy for general structural dynamics problems [27]. Tang and Lou [15] presented a new type of structure-dependent explicit integration algorithm for RTHS using a pole-mapping rule from the discrete domain and named it real-time substructure testing (RST) algorithm. For the RST method, the structure-dependent coefficients are introduced only in displacement increment while those for velocity increment are constant. Kim and Lee [28] also proposed an explicit integration algorithm with controllable numerical dissipation in the high frequency range, and an extended stability limit and improved spectral characteristics are achieved as well. This method has second-order accuracy but is conditionally stable.

Though the performance of the RST method in solving dynamic problems for linear systems has been well studied, its application to nonlinear systems need to be further discussed and evaluated. In this work, the numerical properties of the RST method in RTHS are discussed for nonlinear structural dynamic problems and compared with other two structure-dependent integration algorithms (i.e., CEM and CRM) since they have quite similar formulations and numerical properties. Firstly, the stability and the upper stability limit of the RST method for nonlinear systems with softening and hardening stiffness are discussed using the spectra radius method. Then, the numerical dispersion and energy dissipation characteristics caused by the numerical damping are evaluated. In addition, the overshooting behaviour [25, 29] of the RST method in displacement is addressed as well. Three different simplified engineering structural systems are presented as numerical examples to demonstrate the computational stability and accuracy of the RST method when compared with other algorithms.

2. Formulations of Structure-Dependent Integration Algorithms

2.1. Single-Degree-of-Freedom (SDOF) Systems

The differential equation of motion for a SDOF system under external excitations can be expressed in time domain as [30, 31]where m, c, and k are the mass, viscous damping, and stiffness, respectively; , , and are the acceleration, velocity, and displacement, respectively; is the restoring force of the numerical substructure, where is for linear systems; is the reacting force measured from the experimental substructure; subscript N and E are the numerical substructure and the experimental substructure; and f (t) is the external exciting force. When discretized by time steps, equation (1) can be written in its discrete form:where , , and are the acceleration, velocity, and displacement at ith time step, respectively; is the reacting force of the experimental substructure at ith time step; and fi is the external exciting force at ith time step.

The general expressions of the displacement and velocity increments for the RST, CEM, and CRM methods is found to bewhere ∆t is the chosen integration time step. The coefficients α1, α2, α3, and α4 for the three methods arewhere viscous damping ratio ξ = ξN+ξE, ξN = cN/2n, ξE = cE/2n, , and ωn is the natural frequency of the system determined from the initial system.

Obviously, the coefficients in each method are not all constant, and some of them depend on the initial structural properties (ξ and ωn) and time step ∆t. Substituting equation (6) into equation (5), it is seen that CEM is implicit as there is in the velocity formulation, which is unknown at (i+1)th time step. Substituting equations (4)–(6) into equation (3) giveswhere c and k are the viscous damping coefficient and stiffness of the global structural system, which are equal to (cN + cE) and (kN + kE), respectively. The velocity of the CEM method at (i + 1)th time step can be rewritten in an explicit form as follows:

2.2. Multidegree-of-Freedom (MDOF) Systems

For a MDOF, the coupled differential equation of motion should be written in a matrix form as follows:where M, C, and K are mass matrix, viscous damping coefficients matrix, and stiffness matrix, respectively; , , and are the acceleration vector, velocity vector, and displacement vector at ith time step, respectively; is the restoring force matrix of the numerical substructure, where for linear systems; and and are the reacting force matrix of the experimental substructure and the external exciting force vector, respectively.

For simplicity, C is supposed to be proportional and the orthogonality of modes is available in the following equations. For a MDOF linear system, rewrite equations (4), (5), and (11) in the modal coordinate system as follows:in which, is the modal matrix derived from solving the eigenvalue problem; n is the number of modes; , , and are modal mass matrix, modal diagonal damping matrix, and modal diagonal stiffness matrix, respectively; , , and are the acceleration vector, velocity vector, and displacement vector in modal coordinates, respectively; (j = 1, 2, 3, 4) are the diagonal integration coefficient matrices, and αj is expressed as follows:where and .

For a MDOF nonlinear system, equations (12)–(14) are not valid, and the general formulation and coefficients for the three methods should bewhere ; ; C0 is the damping matrix generally determined from the initial structural properties; and K0 is the initial stiffness matrix.

2.3. Stability Analysis of the RST Method in Nonlinear Systems

The unconditional stability of the RST method for linear elastic systems has been verified by authors [15, 24, 32]. In this section, stability analysis of the RST method for nonlinear systems is mainly carried out and compared with that of the CEM and CRM methods.

In order to realistically reflect the change of stiffness during the test, a parameter, defined as the ratios of the stiffness at the end of each time step ki over the initial stiffness k0, is introduced to monitor this change. That is,where δi+1 is the instantaneous degree of nonlinearity. When δi+1 = 1, it means that the system is linear, and no stiffness change occurs during the test. When δi+1>1, it represents that the instantaneous stiffness becomes harder at the end of the (i + 1)th time step, while 0 < δi+1 < 1 denotes instantaneous stiffness soften.

From equations (3)–(5), the computing sequence at (i + 1)th time step can be written in a recursive matrix form as follows:in which, Xi+1 = [ui+1], is state vector values at the (i + 1)th time step; load factor ; and amplification matrix A for each algorithm is found to be

Based on the algorithm stability analysis theory, a stable computation of an integration algorithm can be obtained when the spectral radius , where the eigenvalues can be determined by solving the following eigenvalue problem [31, 32]:where I is the identity matrix and A1, A2, and A3 are three coefficients indicating half of the trace, sum of principal minors, and determinant of A, respectively.

After substituting equations (21)–(23) into equation (24), the upper stability limits [Ω] for the RST, CEM, and CRM methods can be formulized. The variations of [Ω] for all three methods with δ and the case of ξ = 0.0, 0.05, 0.2 are shown in Figure 1. It illustrates that the RST method is conditionally stable for nonlinear systems except for instantaneous stiffness softening systems when ξ = 0.0. The upper stability limit [Ω] of the RST method increases with an increasing value of δ when 0 <δ < 1 while decreases with a growth of δ value when δ > 1. The initial damping ratio has a great impact on [Ω] as well especially for instantaneous stiffness softening systems, and a larger value of ξ will significantly decrease the upper stability limit of the RST method. The variations of [Ω] with δ for the CEM and CRM methods are also plotted in Figure 1, and both these methods are unconditionally stable when 0 <δ ≤ 1 while conditionally stable when δ > 1. In particular, the curves for the CRM method are overlapped together for different ξ while that for the RST and CEM methods seem to move upward as the value of ξ increases.

3. Accuracy Analysis

Two principal eigenvalues at the (i + 1)th time step for the three structure-dependent integration algorithms are complex conjugate and can be expressed as [33]wherein which and are the numerical damping ratio and frequency related to integration algorithms.

In general, numerical dispersion and energy dissipation characteristics are two indexes evaluating the accuracy of an integration algorithm. The former characteristic is usually expressed by the relative period error , where and referred to the numerical and exact periods of the systems. The later one can be expressed by the numerical damping ratio , which is related to amplitude decay in one cycle of free-vibration.

The variations of PE are shown in Figure 2 (ξ = 0) and Figure 3 (ξ = 0.1) for different δ. It is worth noting that values of PE in all the figures are positive and increase with increasing values of Ω for a given δ, which results in period elongation for the three methods. It is also seen that the curves are totally overlapped for the RST, CEM, and CRM methods for different values of δ as ξ = 0 as well as for δ = 1 as ξ = 0.1. Meanwhile, for a specific value of Ω, the RST method possesses the minimum period elongation while the CRM method produces the maximum period elongation for δ = 0.5, 2 as ξ = 0.1.

The variations of with Ω for the RST, CRM, and CEM methods can be found in equation (26) with A1 and A2 derived from equations (21)–(24). It is noteworthy that for ξ = 0, A2 equals to 1 for all the three methods, which leads to ρ = 1 and  = 0 in equation (26). This implies that no energy dissipation occurred for the three methods if a zero viscous damping ratio adopted. On the contrary, there will be no amplitude decay for an undamped free vibration. When viscous damping ratio is considered, ρ = 1 and the variations of with Ω for different δ as ξ = 0.1 are shown in Figure 4. For a linear elastic system with viscous damping ratio (δ = 1),  = ξ for Ω = 0 and then reduces below ξ with increasing Ω, indicating that the viscous damping ratio is not effective in reducing spurious participation of higher modes. Also, the curves for the three methods are totally overlapped in this case, as shown in Figure 4(b). Meanwhile, drastically different phenomena are found for the case of δ = 0.5 and 2.0. In Figure 4(a), the variations of with Ω are as similar as the case of δ = 1, but  > ξ for a small value of produces undesired dissipation in lower modes. In Figure 4(c), reduces far below ξ with increasing omega and then increases again in a specific value of Ω. It implies that favourable energy dissipation can be obtained to filter out the spurious participation of higher modes.

4. Overshooting

In general, two primary factors are taken into account for choosing an appropriate time step Δt in the step-by-step integration algorithm in order to perform well in solving structural dynamic problems. The first factor is the highest frequency fh that is noteworthy in the dynamic analysis, and its value depends both on the higher modal frequency of the system and the effective high frequency of the dynamic load [15]. The second factor is the stability requirement of an integration algorithm. Since the value of a conditionally stable algorithm is strictly limited by upper stability limit, the value of Δt can be easily determined. If an algorithm is unconditionally stable, the value of Δt depends largely on the first factor.

According to the basic theory of structural dynamics, reliable solutions can be achieved when Δt/T ≤ 1/10 (T = 1/fh) [31]. That is, a large time step involving in dynamic analysis, properly speaking, a large value of the ratio of Δt to T will lead to calculation distortion or so called “overshooting.” Actually, the phenomenon that an integration algorithm might overshoot the exact solution in the initial stage of the transient response or in the steady-state response when Ω ⟶ ∞ has been found and discussed in some papers [25, 29]. Tang and Lou [34] provided a detailed explanation for the issue of overshooting by solving free and forced vibration problems of three single-degree-of-freedom systems with different initial stiffness (k0 = 102, 104, 106) and different values of Δt/T. It is reasonable to discuss the overshooting issue in zero damping systems as undesirable contaminated responses will be eliminated due to vibration damping. It is noted that overshooting in the steady-state response depends mainly on the rationality of integration time step selection. In a real-time hybrid testing, overshooting will not appear in the entire response for any unconditionally stable integration algorithm as the value of time step Δt is strictly limited by the “real-time” requirement. Therefore, the overshooting property of the three methods will not be discussed in detail herein, and the conclusions made by Tang and Lou [34] will be demonstrated by the following numerical examples.

5. Numerical Simulation of RTHS

In this section, three numerical examples are applied to examine the numerical properties of the RST method when compared with the CRM and CEM methods for both linear elastic and nonlinear systems.

Example 1. An undamped 4-story building
A linear 4-story shear building with no structural damping is shown in Figure 5. The initial structural properties are m1 = 108 kg, m2 = 105 kg, m3 = 105 kg, m4 = 103 kg, and k = 360 × 107 N/m. The natural frequencies of the system are 5.9940, 116.95, 306.59, and 1906.9 rad/s for the four modes, respectively; the normalized modal matrix is given by Φ = [ϕ1ϕ2ϕ3ϕ4], where ϕ1 = [0.9970 0.9990 1.0000 1.0000]T, ϕ2 = [−0.0002 0.6293 1.0000 1.0000]T, ϕ3 = [0.0006–1.5882 0.9412 1.0000]T, and ϕ4 = [−0.0000 0.0001–0.0101 1.0000]T. It is clear that the system is intentionally designed to have a relatively high frequency for the forth mode. The first two stories are referred as the experimental substructures while the upper two stories are taken as the numerical substructures.
The free vibration responses of the system subjected to zero initial velocity and displacement U0 = (ϕ1 + ϕ2 + 0.5 ϕ4) are calculated. The time step of Δt = 0.01s is employed for RST, CRM, and CEM, and the ratio of Δt/T4 = 3.036. The displacement response obtained from the undamped displacement equation based on the structural dynamics theory [31] is considered as the exact solution. Defining the relative error e (%) asin which A is the peak value of responses derived from each integration algorithm; subscript “METHOD” refers to the RST, CRM, and CEM methods, while “EXACT” refers to the exact solutions, which are calculated by the Newmark method with a time step of 0.005 s for the overall structural system.
It is shown in Figure 6 that all the three methods provide reliable solutions, and the relative errors of the displacement peaks are 1.46% (RST), 2.64% (CRM), and 0.37% (CEM), respectively. The three methods result in a slight period elongation and zero amplification decay, as shown in Figure 6(a), which are consistent with the previous conclusions. The results also indicate that overshooting behaviours will not appear in a zero damping system if the time step is determined from the perspective of the highest frequency.

Example 2. A 3-story building with a vibration isolator
A linear 3-story shear building with a vibration isolator is shown in Figure 7. In this case, the frame is taken as the numerical substructure while the vibration isolator is supposed to be the experimental substructure. The initial structural properties are m = 160 × 103 kg, kN = 360 × 107 N/m, kE = 9 × 106 N/m, cN = 0, and cE = 115 × 105 N s/m. The natural frequencies of the whole system are 3.7459, 114.10, 212.17, and 277.17 rad/s for the four modes, respectively; the normalized modal matrix is given by Φ = [ϕ1ϕ2ϕ3ϕ4], where ϕ1 = [0.9231 0.9231 1.0000 1.0000]T, ϕ2 = [−1.0000–0.5225 0.4245 1.0000]T, ϕ3 = [1.0833–1.0000 −1.0833 1.0000]T, and ϕ4 = [−1.0000 2.3657–2.3657 1.0000]T.
The seismic performance of the system is studied by exciting the ground acceleration record of EL Centro (1940 NS) with peak ground acceleration (PGA) scaled to 0.8 g. The system is subjected to zero initial velocity and initial displacement. The time step of Δt = 0.02 s is employed for the three methods, and the ratio of Δt/T4 = 0.883. The responses obtained from the Newmark method with constant average acceleration (N4 algorithm) [15] with a time step of Δt = 0.005 s is considered as the exact solution for comparisons. Numerical results for the top story are plotted in Figure 8, and the relative errors of the displacement peaks are 0.37% (RST), 0.50% (CRM), and 0.32% (CEM), respectively.

Example 3. A 5-story building with softening stiffness
A 5-story shear building (as shown in Figure 9) with softening stiffness is applied to investigate the performance of the RST, CRM, and CEM methods in a nonlinear system. In this case, the top three stories are taken as the numerical substructure while the first two stories are supposed to be the experimental substructure. The initial structural properties are m = 105 kg, k0 = 108 N/m, and ξ = 0.05. The softening stiffness of the system after deformation is intentionally designated as , where subscript i refers to the ith story and Δui is the interstory drift for the ith story. The natural initial frequencies of the whole system are 9.001, 26.273, 41.417, 53.206, and 60.684 rad/s for the five modes, respectively; the normalized modal matrix is given by Φ = [ϕ1ϕ2ϕ3ϕ4ϕ5], where ϕ1 = [0.2828 0.5263 0.7368 0.8947 1.0000]T, ϕ2 = [−0.8235–1.1176 −0.5882 0.3160 1.0000]T, ϕ3 = [1.3571 0.3837–1.2143 −0.7143 1.0000]T, ϕ4 = [1.7000–1.4000 −0.5372 1.9000–1.0000]T, and ϕ5 = [−1.8614 3.1643–3.5366 2.6059–1.0000]T. Assuming Rayleigh damping [31], the damping matrix of the system and the proportionality factors are given byThe seismic responses of the system with zero initial conditions is studied by exciting the ground acceleration record of EL Centro (1940 NS) with peak ground acceleration (PGA) scaled to 0.9 g. The results obtained from N4 algorithm with a time step Δt = 0.005 s is considered as the exact solution for comparisons. The time step of Δt = 0.02 s is employed for the RST, CRM, and CEM methods, and the ratio of Δt/T5 = 0.1933. Seismic responses for the 5th story and the shear force-displacement curves for the 1st story are plotted in Figures 10 and 11. It is seen that all the methods can result in acceptable solutions for the nonlinear variation of the stiffness and overshooting behaviours do not appear as well. The relative errors of the displacement peaks for the RST, CRM, and CEM methods are 3.98%, 4.05%, and 3.97%, respectively.

6. Conclusions

Applications of the RST method to nonlinear systems in RTHS have been evaluated and compared with other two structure-dependent explicit integration algorithms, i.e., the CRM and CEM methods, since they have quite similar formulations and numerical properties. Stability analysis indicates that the RST method is only conditionally stable when applied to nonlinear systems, and the instantaneous degree of nonlinearity will exert great influence on the upper stability limit of the method. For stiffness-softening systems, the upper stability limit increases with an increasing value of the instantaneous degree of nonlinearity, while for stiffness-hardening systems, it decreases when the instantaneous degree of nonlinearity becomes larger. Meanwhile, the initial damping ratio of the system has a negative impact on the upper stability limit especially for instantaneous stiffness softening systems, and a larger value of the damping ratio will significantly decrease the upper stability limit of the RST method.

Three simplified engineering structural systems are presented as numerical examples to demonstrate the computational stability, accuracy, and overshooting behaviours of the three methods. Numerical results demonstrate that overshooting behaviours will not appear if the time step is reasonably chosen for integration. It is also illustrated that the RST and CEM method have a better accuracy than the CRM method in solving dynamic problems for both linear and nonlinear systems. However, as the velocity of the CEM method is in an implicit form, the RST method provides more benefits in computational efficiency as both the displacement and velocity are in an explicit form.

Data Availability

The data used to support the findings of this study are included within the article. The results of this paper can be verified according to the data and methods given in the paper.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This paper was based upon work supported by the National Natural Foundation of China (Grant nos. 41904095, 51979027, and 51908048), Natural Science Foundation of Hebei Province (Grant no. E2019210350), Natural Science Foundation of Shaanxi Province (Grant nos. 2019JQ-021), State Key Laboratory of Mechanical Behavior and System Safety of Traffic Engineering Structures at Shijiazhuang Tiedao University (Grant nos. ZZ2020-04), and Fundamental Research Funds for the Central Universities (Grant nos. DUT19JC23 and DUT19RC(4)020).