Rayleigh damping is often used owing to its ease of use. It is a classical damping model expressed by damping terms proportional to mass m and initial stiffness k Eq. 1. c = a 0 m + a 1 k (1) Hence, if coefficients a 0 and a 1 are known, the i th mode damping ratio ξ i can be found using the following expression. ξ i = 1 2 ω i a 0 + ω i 2 a 1 (2) Coefficients a 0 and a 1 can be determined from the specified damping ratios ξ p and ξ q for the p th and q th modes, respectively. For these two modes in matrix form, Eq. 2 becomes. 1 2 [ 1 / ω p ω p 1 / ω q ω q ] { ξ p ξ q } = { a 0 a 1 } (3) These two algebraic equations can be solved to determine coefficients a 0 and a 1 .

The stiffness-proportional term a 1 k of Eq. 1 is a model created on the basis of an assumed initial stiffness regardless the nonlinearity of response history analysis (RHA). However, Chrisp (1980) reported that an unintended spurious damping force is generated in plastic hinges, and the mechanism of its generation was explained later by Bernal (1994). This phenomenon occurs when a stiffer nonlinear element is explicitly incorporated into the end of a beam. Additionally, Léger and Dussault (1992) and Charney (2008) pointed out nonlinearity may change the natural frequency. Therefore, to avoid all aforementioned problems, a model in which the stiffness term is proportional to the tangential stiffness k T rather than the initial stiffness k used. c ∗ = a 0 m + a 1 k T (4)

The damping force f d is calculated as f d = c ∗ ⋅ x ˙ using from the tangential damping coefficient c ∗ and velocity x ˙ . Assuming the spring is completely bilinear, its behavior exhibits initial stiffness proportionality in the elastic stage and drastically the damping force suddenly becomes zero after yielding. This behavior lacks a physical basis and this discontinuity in the damping force poses problems in numerical calculations.

Wilson and Penzien (1972) expresses the modal damping matrix C W P by the following equation. C WP = M ( ∑ q = 1 n 2 ξ q ω q M q φ q φ q T ) M (5) Where, M is the mass matrix; and ξ is the damping ratio; φ q and ω q are the undamped natural vibration mode vector and natural frequency of the q th mode, respectively; and M q is the generalized mass φ q T M φ q . From Eq. 5, if the rotational inertia mass is zero, no damping force is generated in the rotational degree of freedom. Therefore, the spurious damping force generated by the plasticization of the beam-end-moment is theoretically not generated. Chopra and McKenna (2016) suggested using constant modal damping as a countermeasure to the spurious damping force, but it leads to a high computational load because it uses a dense matrix. To avoid this problem, Chopra and McKenna (2016) proposed using tangential stiffness in solving simultaneous equations and treat the difference from the correct mode damping force as an unbalanced force in the equation of motion. Then, this idea was incorporated into OpenSees (McKenna, 1997). However, similar to the problem of Rayleigh damping discussed by, e.g., Hall (2006) and Ryan and Polanco (2008), Pant and Wijeyewickrema (2012), Pant et al. (2013), Hamidreza et al. (2019), if the foundation is not fixed, as in the case of an uplifting structure or base isolation system, or if there is a large peeling between elements such as the formation of cracks in concrete, the deformation will be underestimated, and the model may become nonconservative.

In addition, Luco and Lanzi (2019) reported that in modal damping, after the nonlinear element transition from the elastic state to plastic state, the degree of freedom without mass becomes an unintended velocity response because of numerical artifacts and the absence of damping terms. To avoid this problem, model damping can be improved by adding infinitesimal stiffness-proportional damping.

The capped viscous damping model is a damping model that can be applied to a complex vibration model with large nonlinear behavior; for example, it can be applied to the sliding or uplifting of a foundation. Although tangent Rayleigh damping is effective in such an analysis, there is concern that tangent Rayleigh damping has little physical basis for reducing the damping force and that damping is overestimated owing to the influence of the mass term. Hall (2006) discusses the increase in the ratio of damping force to spring strength using a 10-story lumped mass shear model as an example. In the elastic state, the ratio of the restoring force and damping force is kept constant, but after yielding, the restoring force is limited, whereas the damping force increases in proportion to the velocity. This is unnatural and cannot be explained as a real phenomenon. In the case of initial-stiffness-proportional damping, the restoring force f k and damping force f d due to the arbitrary story stiffness k have the following relationship. f k + f d = ( 1 + 2 ξ i ) k x (6)

If the spring is elastic and oscillated in the first mode, the peak values of the restoring force and damping force will maintain a ratio of 2 ξ with a phase difference of π/2. However, the spring exceeds the yield displacement and the restoring-force peaks at yield strength R m a x . Hence, the ratio to the damping force exceeds 2 ξ . Hall (2006) reports that the computed maximum value of damping force was 60% of the yield strength of the building, and this response result is nonconservative. Therefore, this problem can be avoided by setting the capped damping force as shown below [Hall (2006)]. f d = ± min ( a k k | x ˙ | 、 2 ξ R m a x ) (7)

Here, R m a x should be set before the analysis is executed, and there is no established idea regarding this at present. Therefore, it is necessary to set R m a x individually depending on the designer’s engineering judgment. In previous works Hall (2006), Hall (2018), and Qian et al. (2021), as examples of setting R m a x , the yield strength of the story is assumed and determined by calibration.

In addition to the effect of suppressing an abnormal increase in the damping force after yielding, caped damping has the characteristic of causing a certain degree of frequency insensitiveness. The capping fraction 2 ξ is based on the ratio of the damping and stiffness forces when a linear structure is vibrating in a mode at the resonant frequency of that mode. The single-degree-of-freedom oscillator with capped viscous damping is illustrated Figure 1. If the system is oscillating at a given displacement amplitude ( X a ) with a random frequency below the value of ξ , it will behave as an initial-stiffness-proportional damping oscillator (red line). Further, if the system is oscillating at the same X a but with higher frequencies, it will be characterized with more damping than desired. However, a higher-mode damping does not increase with increasing frequencies due to a capping effect. Currently, the damping force of the oscillator has approached a frequency-insensitive state. Therefore, the capped damping force is largely dependent on the frequency-insensitive behavior. To further enhance the frequency-insensitive nature, we evaluated an adjusted capped value according to the restored force amplitude amount, which was obtained from moment to moment in some way. This approach was better than fixing the capped value in advance. Furthermore, according to previous studies, the analyzed targets have been limited to a one-way input; therefore, it is necessary to model the system so that it can be applied to a three-way simultaneous input.

The limitations of the equivalent viscous damping model identified in previous studies are summarized as follows:

1) Generally, the viscous damping model depends on frequency (Clough and Penzien (2003)).

In this study, we focused on the capped viscous damping proposed by Hall (2006) with the aim of exploring a well-balanced model that can avoid the aforementioned problems. The vibration characteristics were analyzed by comparing it with several damping models, and the engineering convenience/effectiveness of this damping model was considered.

If the building undergoes simple vibration in the first mode and the displacement is the maximum amplitude X a , the maximum damping force D m a x can be expressed using the story stiffness k. D m a x = 2 ξ ⋅ k ⋅ X a = 2 ξ ⋅ R m a x (8) R m a x is the maximum restoring force, which is equal to k ⋅ X a . And the building is oscillated in the first mode, the damping will be linear viscous damping unless the restoring force is capped by the yield strength. Next, considering the vibration system oscillates simply in the second mode, the maximum damping force can be similarly evaluated by the following equation using the maximum amplitude X a . D m a x = 2 ξ ω 1 ⋅ k ⋅ ω 2 ⋅ X a = 2 ξ ⋅ ω 2 ω 1 ⋅ R m a x (9)

Equation 9 shows that the damping force increases in proportion to ω 2 / ω 1 , Where ω 2 is the second mode natural frequency, and this is a characteristic of the initial-stiffness-proportional damping. However, the damping force is capped at 2 ξ R m a x , the system vibrating in the second mode has a frequency-dependent coefficient ω 2 / ω 1 of 1.0. Considering that the second mode vibrates with the same amplitude, the amount of damping energy dissipated per cycle is about the same. If this principle can be incorporated into numerical analysis well, it will be possible to realize modeling that does not depend on the frequency but depends on the amplitude. Therefore, it is necessary to determine the damping force that depends on the amplitude of the current restoring force. However, a phase difference of π/2 occurs at the time of the maximum response of displacement (i.e., restoring force) and velocity. As a simple method, it is conceivable to determine the capped damping force with the maximum restoring force experienced in the past as R m a x , but the frequency insensitiveness performance may be inferior after the maximum value experience. R m a x to be set during the time history analysis is an important topic. Therefore, in this study, it is assumed that the relationship of X ˙ a ≅ ω X a is approximately established between the maximum displacement X a and velocity X A ˙ recorded near an arbitrary time, and the capped damping force D m a x can be evaluated by D m a x = 2 ξ ⋅ R m a x using the maximum restoring force R m a x recorded most recently. This is based on the rule of thumb that the maximum displacement X a and velocity X A ˙ near a certain time are often governed by one primary mode with a natural frequency ω . Therefore, the secondary higher-mode damping force cannot behave independently because of the capped damping force D m a x . Hence, it is considered that frequency insensitiveness can be achieved to some extent. However, this concept may have weak points in that the prediction accuracy of the maximum damping force value is poor in the case where multiple modes are overlapping, e.g., in the case of simultaneous input RHA in three directions and RHA of high-rise buildings where higher modes are prominent. This point is discussed in the Discussion section.

In the next section, we explain the setting of the maximum restoring force R m a x .

This study proposes a capped damping force 2 ξ R m a x ; here, R m a x is the most recently recorded restoring-force peak. Figure 2 shows the time history of any restoring force at the nodal end of the stiffness element. The procedure for updating R m a x will be explained using this restoring-force time history as an example. R m a x stores the maximum value in each of the positive and negative directions and updates it according to the following rules.

1) If R m a x is compared with the current time restoring force and if the current time restoring force is larger, R m a x is updated by the current restoring force.

The red and blue dashed lines in Figure 2 indicate the time domain in which R m a x on the positive and negative sides is updated, and circle indicates the time when R m a x was updated at the end immediately before the decrease. In the figure, “A” represents the state where the restoring force is increasing to the positive side, and at this time, the peak value R p m a x on the positive side is continuously updated along the red dashed line according to Rule I. “B” represents the stage of switching to the negative side, and the restoring force increased to the negative side, and the peak value R n m a x on the negative side was sequentially updated by Rule I in the area of this blue broken line. “C” represents the stage of increase on the positive side again, and R p m a x is sequentially updated by Rule I in the area of the red dashed line that exceeds the most recently recorded R p m a x . “D” represents the stage of increase on the negative side again, and although it does not exceed the most recently recorded R n m a x , the value of R m m a x is updated by Rule II. Even if the direction of the restoring-force changes from decreasing to increasing on the positive side as at “D’,” the R n m a x is not updated.

As described earlier, by setting the peak value of the most recently recorded restoring force to R m a x , a damping force that depends not on the frequency but on the amplitude near the current time is created. There is no distinction between the elastic and inelastic elements in the update of R m a x , and the restoring force at that time is applied to the judgment. By updating the R m a x of each restoring-force vector of all elements according to the above rules, the general-purpose application in three-dimensional dynamic analysis will become possible.

Figure 3 shows the program flow which is written specifically test capped viscous damping. For the total equation of motion, incremental Eq. 10 at the integral time interval Δ T is cumulatively obtained. M Δ x ¨ + C Δ x ˙ + K Δ x = Δ p (10)

Here, M is mass, K is tangent stiffness, and C is the initial-stiffness-proportional damping matrix. The incremental equation forms the form of the relational expression Eq. 11 between the static load and displacement. The effective stiffness matrix K ˜ is as shown in Eq. 12, and the effective load Δ p ˜   is as shown in Eq. 13. K ˜ Δ x = Δ p ˜ (11) K ˜ = K + 6 Δ T 2 M + 3 Δ T C (12) Δ p ˜ = Δ p + ( 6 Δ T x + 3 x ¨ ) M + ( 3 x ˙ + Δ T 2 x ¨ ) C (13)

The program performs the triangular decomposition of the effective stiffness matrix K ˜ by the modified Cholesky decomposition in 1) as shown in Figure 3. Next, the effective load shown in Eq. 13 is calculated by 2). The element stress and maximum restoring force R m a x are updated in 3). Finally, in 4), the smaller of the element damping force vector f d obtained by the element initial-stiffness-proportional damping matrix c in Eq. 14 and the element damping force vector obtained from the maximum restoring force is adopted as the element damping force f d ' . Further, i in Eq. 15 represents the degree of freedom of the local element. Here R m a x is a two-dimensional array and holds the maximum values of the positive and negative sides of the restoring force to evaluate the positive and negative directions of the f d , respectively. f d = c ⋅ x ˙ ,     c = 2 ξ ω 1 ⋅ k (14) f d ′ ( i ) = min ( f d ( i ) ,   2 ξ R m a x ( i ) ) (15)

The balance of the total equation of motion is calculated using the updated damping force vector d ' which is the sum of each local element damping force vector f d ' and the unbalanced force u is calculated using Eq. 16 and carried over to the next step. p and r are the external force and restoring-force vectors at the current time. u = p − M x ¨ − d ′ − r (16)

To evaluate the differences in the responses of different damping models, a comparative study was conducted using a fish bone model that simulates a high-rise building and a period for 2 s (Figure 4A). The vertical degrees of freedom of all nodes were constrained, and the horizontal displacement of the nodes at each beam level were equal to the that of the beam–column joint node by multi-point constraints. The story height was 3.5 m, and the span is assumed to be 6 m, and 3-m-long beams were set on the left and right. The beam and column sections are listed in Table 1. The beam stiffness was multiplied by Φ (=1.5) considering the slab, and the column stiffness was about 1.8–1.9 times that of the beam. A mass (59.19 t) was uniformly applied to each floor at the beam–column joint node as shown in Figure 4A so that the period was 2.0 s and the rotational inertia mass of each node was zero. To confirm the differences depending on the types of structures, two cases were considered: the case of the beam-end assuming a steel (S) structure with restoring-force characteristics is shown in Figure 4B, and the case of the beam-end assuming a reinforced concrete (RC) structure with the Takeda model is shown in Figure 4C. The initial stiffness and yield moment of the RC model are the same as those of the S model for convenience. The crack-bending strength of the Takeda model is assumed to be 0.3 times its yield strength, and the post-cracking stiffness is assumed to be 0.5 times the initial stiffness. Table 2 lists the initial stiffness (3EI/L) of the beam, post-cracking stiffness ratio, cracking moment, and yielding moment. The inelastic spring at the end of beam is generally incorporated into the element stiffness matrix implicitly or arranged explicitly with an additional node. In the arranged explicitly with an additional node case, the initial stiffness of the end hinge spring is set at a higher value so as that the linear natural vibration of the model is unaffected; however, the spurious damping force is prominent in such modeling (Chrisp (1980); Bernal (1994)). In this study, this method is used to insert an inelastic spring with stiffness 1,000 times higher than the beam (Figure 4D). Therefore, for the inelastic spring with nonlinear characteristics (Table 2), K₁ should be multiplied by 1,000, and conversely, the post-cracking stiffness ratio of K₂ and K₃ should be multiplied by 1/1,000, and K₂ gives 0.0005 and K₃ gives 1.0E-6.

Table 3 lists the natural period and effective mass ratio (cumulative values from lower mode).

Ground motion was created by using an acceleration response spectrum with 5% damping as the target spectrum and using the random phase characteristic (return period of 500 years). Figure 5 shows the simulated ground motion and acceleration/pseudovelocity response spectra.

The three types of capped viscous damping models, including the original capped viscous damping [CP(O)], capped viscous damping A [CP(A)], and capped viscous damping B [CP(B)] models, were compared to prove the effectiveness of the proposed modeling. The accuracy of each damping model was verified by comparison with the nonlinear constant modal damping model (WP), which creates a damping matrix by performing an eigenvalue analysis each time the structural stiffness changes.

The original or CP(O) model has been proposed as a method for giving a lateral interstory damper (Hall (2006); Hall (2018)). The same modeling was used in this study. The height-wise distribution of damper coefficients ( C i ) was taken as proportional to the lateral interstory stiffness ( K i ), i.e.,   C i = α K i (where i is a story number). With this assumption, the coefficient ( α ) for the assumed first mode damping ratio ( ξ 1 ) was calculated by Eq. 17 Qian et al. (2021). ξ 1 = α T 1 ∑ i K i ( ϕ 1 , i − ϕ 1 , i − 1 ) 2 4 π ∑ i m i ϕ 1 , i 2 (17) Where, m i , ϕ 1 , i , and T i denote the mass lumped at the floor, lateral displacement vector of the first mode, and first mode period, respectively. K i was obtained by static elastic analysis with a horizontal external force ( ϕ 1 , i ) distribution. Based on these parameters, α = 0.0127 . A story strength ( R i ) value was estimated as the maximum story shear force responses, which resulted from a nonlinear RHA due to the assumed seismic motion with ξ = 2 % of the nonlinear constant modal damping. The results estimated from the measurements of elastic and inelastic RHA in RC and S structures are summarized in Table 4. The damper coefficient after the relief velocity is 1/1,000 times the initial damper coefficient.

The CP(A) is a method for determining the capped damping force. This was done by setting the maximum value of the restoring force experienced in the past as R m a x . This model was not the recommended model in this study; however, it was considered as the target for comparison because of its clearer mechanism. This model was also easier to implement in the program used than CP(B).

The CP(B) was considered as the main method in this study. This is a method for determining the capped damping force from the maximum value of the most recently experienced restoring force, which was set as R m a x .

The oscillation characteristic of constant modal damping WP, CP(A), and CP(B) were compared using three cycles of sinusoidal wave synchronized with the building period; further, CP(O) was not compared. The input wave is shown in Figure 6A. Using this input wave, the steady vibration at the time of input and the free vibration after the input can be examined in one analysis. The maximum acceleration was 100 gal, and the response analysis time was 30 s. The integration step was 0.001 s; the damping ratio was 2%; The responses obtained by using each sinusoidal wave tuned to the first mode to the third-mode periods of the frame were compared. Figure 6 shows the maximum story drift angle and acceleration for each floor. Figure 6B is the result of the first mode periodic input. The maximum responses of CP (A) and CP (B) consistently show the same results in the first to third modes. Further, CP(A) and CP(B) have slightly larger responses than WP. This is because as suggested in the previous section, the velocity precedes the displacement, and hence, the damping force peaked at the maximum value of the past restoration force, and the damping force could not be fully exerted. Moreover, as shown Figures 6C,D, the maximum response of the CP(A) and CP(B) of the second- and third-mode periodic inputs was the same as that of the WP. Figure 7 shows the mode RHA, which is the mode response y s extracted by the expansion theorem shown in Eq. 18. y s = { ϕ s } T [ M ] { x } { ϕ s } T [ M ] { ϕ s } (18) { x } : Response displacement vector { ϕ s } : sth mode eigenvector.

From Figure 7, we can see that the maximum response occurs at the final peak of the sweep excitation, but CP(A) and CP(B) show a relatively frequency insensitiveness characteristics for the response of such an amplification process. However, in the free vibration after the forced-excitation phase, there is a tendency to overestimate the damping slightly compared to the estimates of WP. This overestimation is attributed to the fact that since R m a x is the previous peak value, in the case of free vibration in which the amplitude gradually decreases, the R _max is evaluated to be larger than the amplitude at the current time. CP(B) did not significantly increase the damping during free vibration compared to CP(A) because the R m a x evaluation continued to be updated even after the maximum peak.

Next, Evaluate the oscillation characteristics of real-phase seismic waves as elastic and inelastic. To compare the differences in the characteristics of the steel and reinforced concrete structures, the integration time step was 0.0001 s, and the damping ratio was 2% for both S and RC structures.