Error Caused by Damping Formulating in Multiple Support Excitation Problems

Featured Application: The conclusions acquired in this paper can be used for assessing damping modeling errors in multiple support excitation problems. The modeling method proposed in this paper can be used for both spectral analysis, linear or nonlinear time history analysis of large span structures. Abstract: The e ﬀ ect of multiple support excitation is an important issue in studying large-span structures. Researchers have shown that the damping related terms in the equation of motion can induce errors in the analysis. Wrongly modelling the damping matrix can induce false damping forces between the structure and the reference coordinates. In multiple support excitation problems, this error is increased when absolute coordinates are used. In this paper, this part of the error is deﬁned as virtual damping error. The error caused by using Rayleigh damping instead of Modal damping is called damping truncation error. This study focuses on the virtual damping error and the damping truncation error that exist in the modeling methods widely used in multiple support excitation problems, namely, large mass method (LMM), relative motion method (RMM), and absolute displacement method (ADM). A new Rayleigh damping formula is proposed for LMM to prevent virtual damping error. A form of equation of motion derived from the converged LMM was proposed in the authors’ previous work. This equation of motion is proved in this paper to be equivalent to RMM when modal damping and the new Rayleigh damping formula are used. RMM is proved free from the virtual damping error. The inﬂuence of multiple support excitation e ﬀ ect on the damping formulating errors is studied by spectral analysis. One simpliﬁed spring-mass model and two bridge models are used for numerical simulation. The results from the numerical simulation testify to the conclusions from the spectral analysis.


Introduction
The effect of multiple support excitation on the responses of large-span structures has been studied for decades. Errors that existed in modeling structures have been studied extensively. There are three methods mainly utilized for modeling the structures in multiple support excitation (MSE) problems: Relative motion method [1] (RMM), large mass method [2] (LMM) and absolute displacement method [3] (ADM).
RMM is the most widely used method. In RMM, the response of the structure under multiple support excitation is separated into dynamic responses and pseudo-static responses, these two and is used as a reference to measure the errors that exist in the other methods. This EOM is further extended to the Rayleigh damping method. The error ratios of the resonance response of the modal frequency response function derived from different modeling methods are derived. Next, the influence of MSE on both the virtual damping error and the damping truncation error is discussed with the random vibration method. Lastly, cascaded spring-mass systems of different degrees of freedom, a 4-span rigid frame bridge and a large span cable-stayed bridge are used to study the error induced in the modeling method and demonstrate the influence of MSE on the damping modulating error.

Review of Modeling Methods
The free vibration EOM of a structure containing the excitation degrees of freedom is shown in Equation (1).
The matrices M, C, K are the mass, damping, and stiffness matrices, respectively, for different degrees of freedom; The subscript s and b refer to the structural degrees of freedom and the excitation degrees of freedom; Y stands for the absolute displacement of degrees of freedom and 0 stands for the zero matrices.
In RMM, the responses of a structure are separated into dynamic responses and pseudo-static responses as shown in Equation (2).
The pseudo-static responses satisfy the following equation: Define R = −K −1 ss K sb then the pseudo-static responses can be solved by Equation (4).
The dynamic responses of the structure can be solved by Equation (5).
According to Clough et al. [1], the damping term at the right-hand side of Equation (5) is much smaller than the inertial term and thus is negligible and the dynamic responses of the structure can be solved by Equation (6). Most of the studies considering multiple support excitation adopted Equation (6) as the EOM of the dynamic responses.
The EOM of LMM is shown in Equation (7), where M bb stands for the added large masses.
Appl. Sci. 2020, 10, 8180 4 of 26 Extract the first line of Equation (7) and ignore the inertial and damping contribution of the support. The EOM of ADM is shown in Equation (8).
As has been mentioned in the introduction, the error of the aforementioned methods mainly comes from the damping formulating process.

Error in the Frequency Response Function of Different Modeling Methods
In this section, the virtual damping error and the damping truncation error in the three modeling methods, namely, LMM, RMM, ADM, are studied by inspecting the frequency response function (FRF) of the modal participation factor in the frequency domain.
Abbreviations such as LMM-R and LMM-M are used to represent using the Rayleigh damping method in LMM and using the modal damping method in LMM, respectively; UN is used to represent the modelling method used in uniform excitation problems. In this section, for ease of analysis, the assumption of lumped mass is used, which is widely used in the field of civil engineering.
To increase the readability of this paper, some symbols used in this article are listed as follows. The response is symbolized as y Lm , where the superscripts represent the modeling method used in deriving the structural EOM; the capital letter L stands for the modeling method used and the lowercase m stands for the damping formulating method is modal damping method; the lowercase r in the superscript stands for the damping formulating method is the Rayleigh damping method. The other Capital letters are listed in Table 1. The method used for uniform excitation problem in relative coordinates

Error in LMM
Properly formulating a support damping matrix can eliminate the virtual damping error. Qin et al. [8] derived the support damping matrix from Equation (7) with the modal damping method.
In the modal damping method, the modal damping matrix C ss C sb C bs C bb is calculated with M ss M sb M bs M bb + M bb and K ss K sb K bs K bb ; Qin et al. [8] showed that when the large masses approach infinity, C ss in C ss C sb C bs C bb converges to the modal damping matrix calculated with M ss and K ss , at the same time, the support damping matrix C sb converges to where → stands for convergence. When C sb = C ss K −1 ss K sb , Appl. Sci. 2020, 10, 8180 5 of 26 where Φ sb T Φ bb T T stand for the rigid modes and the modes related to the relative movement of the supports. The damping matrix in Equation (7) satisfies Equation (10), which is to say, in LMM when the modal damping method is used and when the large masses approaches infinity the virtual damping equal to zero.
C sb can be seen as a damping compensator that can eliminate virtual damping from the modal of the upper structure. It can work with the other support damping in soil-structure interaction problems.
Qin et al. [8] introduce an integrated EOM as shown in Equation (11) by extracting the first line of Equation (7) and assuming C sb = C ss K −1 ss K sb .
This EOM is equivalent to LMM when the large masses approach infinity and has a simpler expression compared with RMM. In the next subsection, this equation is proved to be equivalent to RMM.
In the following context of this section, the EOM derived by the LMM-M would be used as a reference to show the damping truncation error in the LMM-R method.
In linear systems, the responses of a structure can be expressed as the linear combinations of all the modal responses. The damping truncation error will be studied in the frequency domain of modal coordinates. The modes of the structure in Equation (7) is divided into two parts. The first part contains the modes having zero support motion, called structural modes in this paper; the others are called excitation modes. Transform Equation (7) into modal coordinates, the dynamic equation for each modal participation factor is shown in Equation (12). The mode shape matrix Φ ss is normalized to satisfy Φ T ss K ss Φ ss = 2I, where I represents the identity matrix. The process of derivation is shown in Appendix A. ..
where y Lm i is the response of mode i; ω i and ζ i is the natural frequency and damping ratio for the i th mode, respectively; Φ ss T K sb i, * stands for the i th row of the matrix Φ ss T K sb . Equation (12) is actually the same as the dynamic part of Equation (6) in modal coordinates. The responses of the excitation modes are the counterpart of the pseudo-static part of RMM in LMM; and this part can be solved directly from Equation (4) according to Equation (A6) and the conclusion that Φ bs → 0 . This derivation shows in one aspect that LMM is equivalent to RMM when the large masses approach infinity. This conclusion also justifies the assumption of the pseudo-static motion in RMM. Therefore, only the structural modes are studied. When the Rayleigh damping method is used in LMM, the damping matrix should be calculated without the large masses, as shown in Equation (13), otherwise, the large mass proportional damping would affect the accuracy of the support motions.
α 0 and α 1 are parameters of the Rayleigh damping method. In the Rayleigh damping method, the damping ratios of the two modes with frequency ω r 1 and ω r 2 are determined at first. α 0 and α 1 are then calculated to make sure the damping ratio of these two modes being the expected value ζ r 1 and ζ r 2 . When these two damping ratios are assumed to be equal, i.e., ζ r 1 = ζ r 2 = ζ, α 0 and α 1 can be solved by Equation (14).
The damping ratios of the other modes are α 0 /2ω i +α 1 ω i /2. The first mode of the structure is usually chosen as ω r 1 ; the modes whose natural frequency larger than ω r 2 would have a damping ratio Appl. Sci. 2020, 10, 8180 6 of 26 larger than ζ so that the energy of these modes is underestimated. Therefore, ω r 2 is chosen carefully that the underestimation of these modes would not significantly affect the result. For this reason, in this paper, only the damping truncation error of the modes whose natural frequencies are between ω r 1 and ω r 2 are analyzed. The relationship among ω r 1 , ω r 2 and the natural frequency of some of the modes are briefly depicted in Figure 1.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 27 The damping ratios of the other modes are α ω α ω The first mode of the structure is usually chosen as ω 1 r ; the modes whose natural frequency larger than ω 2 r would have a damping ratio larger than ζ so that the energy of these modes is underestimated. Therefore, ω 2 r is chosen carefully that the underestimation of these modes would not significantly affect the result. For this reason, in this paper, only the damping truncation error of the modes whose natural frequencies are between ω 1 r and ω 2 r are analyzed. The relationship among ω 1 r , ω 2 r and the natural frequency of some of the modes are briefly depicted in Figure 1.
To compare Equation (12) The error i z can be expressed in the frequency domain as follows.  The smallest possible value is min ζ r i = 2ζ √ ω r 1 ω r 2 /(ω r 1 + ω r 2 ), when ω i = √ ω r 1 ω r 2 and ζ r i is the damping ratio of the i th mode. Usually, the smallest value cannot be achieved.
Transforming the first row of Equation (7) into modal coordinates, the EOM is shown in Equation (15).
To compare Equation (12) with Equation (15), a variable z i = y Lr i − y Lm i is defined to measure the error when the Rayleigh damping method is used in LMM. Defining another parameter, The error z i can be expressed in the frequency domain as follows.
where z i stand for the Fourier transform of z i . Define µ i = ω i /ω r 1 and ν = ω r 2 /ω r 1 . In the vicinity of zero frequency, where is a small positive number. In Equation (18), the error transfer function at zero frequency is zero; the amplitude of this error increase with frequency ω, and is negatively correlated to the fundamental frequency of the structure, ω r 1 ; the larger the damping ratio ζ is the larger the error is.
Define the error ratio r Lr i = y . The value of r Lr i at zero frequency is infinite, which attribute to the second term at the right-hand side of Equation (5). As shown in Equation (17), the error has a peak in the vicinity of the resonance frequency. The value at the resonance frequency is derived as follows, where η i = ζ/ζ r i , η i represents the multiple of the expected damping ratio over the damping ratio of mode i.
In Equation (19), for a particular structure, the first term in the square root symbol increases with the parameter η i and the second term decrease with µ. The possible largest value of η i is shown in Equation (20), and is achieved in the mode having the smallest damping ratio.
From Equation (19), the behavior of the peak of r Lr i is complex, the natural frequency of the mode having the largest error is between ω r 1 and √ ω r 1 ω r 2 . The largest value of r Lr i is larger than the larger value in r Lr i ω i =ω r 1 = 2ζ and r Lr indicates that, for a particular structure, the larger the value of ν is the larger the error is. The value max(η) doesn't have an upper limit, as a result, the error ratio r Lr i ω i can be a very large value.
When ω approaches infinity r Lr i ω→∞ = 0. It is obvious in Equation (17) that the error ratio would be smaller than r Lr i ω=ω i for the frequencies larger than the resonance frequency. As a consequence, the error ratio at resonance frequency is a preferable value to represent the error in LMM-R. Table 2 shows examples of the value of r Lr for different parameters. In Table 2, r Lr i increases with the increase of ν and ζ; the influence of ν is more significant than ζ. The error ratio can be as large as 43.3%. The value in Table 2 is a good reference for estimating the largest damping truncation error of a particular mode.
An example is demonstrated in Figure 2 showing the shape and the trace of peaks of r Lr i when the mode has different natural frequencies. In Figure 2, the large value of r i L r near 0 Hz is aroused by the second term at the right-hand side of Equation (15). As is analyzed, though the error ratio near 0 Hz is very large, the absolute value of error is small. The large value is the main component of virtual damping error in LMM-R. When ω < ω i , the error ratio decrease to zero at first. After that, the error ratio increase and reach its peak at the In Figure 2, the large value of r Lr i near 0 Hz is aroused by the second term at the right-hand side of Equation (15). As is analyzed, though the error ratio near 0 Hz is very large, the absolute value of error is small. The large value is the main component of virtual damping error in LMM-R. When ω < ω i , the error ratio decrease to zero at first. After that, the error ratio increase and reach its peak at the resonance point, showing that the influence of the large value around 0 Hz and the errors around the resonance point are separate. As analyzed before the peak of the error ratio reach is top when ω i = √ ω r 1 ω r 2 , which is the main feature of the damping truncation error of the Rayleigh damping method. Studies have shown that in the MSE problem, more modes need to be considered [15]. The damping truncation error might not be negligible in the MSE problem.
In LMM-R, the damping matrix C does not satisfy Equation (10). This means that the damping truncation error of LMM-R contains the virtual damping error. At the end of the next subsection, another Rayleigh damping formula is introduced, that can eliminate virtual damping in LMM-R.
In the next subsection, using RMM with the Rayleigh damping method is proved free from virtual damping and has the same form as Equation (11). LMM has its advantage in time history study. The free-from virtual damping version of LMM can be easily derived from Equation (11). In this version, the Rayleigh damping is shown as follows.
where M sb = M ss K ss −1 K sb . With the help of M sb , the second term at the right-hand side of Equation (15) is canceled. The EOM of the modal participation factor of mode i is shown as follows. ..
Except for the damping ratio, the equation has the same form as Equation (12). With Equation (22), the EOM of LMM-R approaches RMM-R when the large mass approach infinity. With the help of Equation (22), the EOM of the LMM-R would have no virtual damping error but has only the damping truncation error inside.

Error in RMM
The way RMM solves the response of the structure guarantees that the EOMs derived with RMM does not contain virtual damping. In Appendix A, the following equation is proved.
Equation (24) indicates that the mode shapes related to the relative motion of the supports are not affected by the mass distribution of the structure. This conclusion gives the physical meaning to separate the responses of the structure into pseudo-static responses and dynamic responses in RMM. Substituting Equation (24) and (4) into Equation (2), Equation (25) can be obtained.
In Equation (25), the dynamic responses equal to subtracting the influence of the rigid body modes and the modes related to the relative motion of the supports from the total responses. In other words, RMM does not contain virtual damping.
To show this property of RMM in another aspect, put all the equations in RMM in a single equation. Adding Equation (4) to Equation (6) and transforming the result into absolute coordinates, the EOM of the absolute displacement is Equation (26). (26) is the same as the one introduced by Qin et al. [8], which is derived from a modified LMM. The form of C sb in Equation (9) guarantees the establishing of Equation (10). The integrated EOM in Equation (26) is more efficient and can also be used in nonlinear analysis.
The above derivation proved that the assumption in RMM that C sb + C ss R = 0 will not induce the virtual damping error regardless of the damping formulating method.
But using the Rayleigh damping method would induce the damping truncation error in the calculation. Transforming Equation (26) into modal coordinates, the dynamic equation for the modal participation vector is shown in Equation (27).
where s i = ω/ω i , y R i is the modal participation factor of mode i and y R i is the Fourier transform of y R i ; ζ i is the damping ratio this mode. the damping truncation error is measured with error ratio r Rr i , the EOM derived with RMM-M is used as the reference. When s i = 1, The superscripts Rm and Rr represent the modelling method RMM-M and RMM-R, respectively. In Equation (28), when ζ is small enough, it has little influence on r Rr i s i =1 ; r Rr i s i =1 increases with η; the largest value of r Rr i s i =1 is obtained in the mode whose natural frequency is the closest to the frequency √ ω r 1 ω r 2 . When s i = 2, When s approaches infinity, The larger η i is, the larger the error is. For a particular mode, the error ratio decreases rapidly from the resonance frequency and then increases gradually as s approach infinity. The increasing rate is very small, r Rr i s i =1 is a good reference to the magnitude of the damping truncation error in RMM-R. Table 3 shows the error ratio for different parameters. It is obvious in Table 3 that the influence of ν is larger than s i and ζ. For a particular structure, the more modes need to be considered the larger the error r Rr i is. Figure 3 is an example used to show the feature of r Rr i and its variation with the change of the modal natural frequency; the parameter settings are the same as Figure  It is obvious in Table 3 that the influence of ν is larger than i s and ζ . For a particular structure, the more modes need to be considered the larger the error r i R r is.  As is shown in Figure 3, r i R r has its first peak at the resonance frequency and the largest peak is acquired when When ω > ω i , r i R r decreases rapidly at first and gradually increases afterward, it finally arrives at a limit as shown in Equation (30). The largest value of the limit is also acquired when ω ω ω = 2 1 i r r . Though the limit error ratio can be very large, the absolute error is small.
In Figure 3, the value of r i R r when ω ω = 2 i are marked with a circle on each line, at these points the value of the modal frequency response function has reduced to 1.9% of its resonance peak; and the value of the modal frequency response function is monotonously decreasing until infinity. As is shown in Figure 3, r Rr i has its first peak at the resonance frequency and the largest peak is acquired when ω i = √ ω r 1 ω r 2 . When ω>ω i , r Rr i decreases rapidly at first and gradually increases afterward, it finally arrives at a limit as shown in Equation (30). The largest value of the limit is also acquired when ω i = √ ω r 1 ω r 2 . Though the limit error ratio can be very large, the absolute error is small.

Error in ADM
In Figure 3, the value of r Rr i when ω = 2ω i are marked with a circle on each line, at these points the value of the modal frequency response function has reduced to 1.9% of its resonance peak; and the value of the modal frequency response function is monotonously decreasing until infinity.

Error in ADM
In this section, the error caused by ignoring support damping, i.e., the ADM, will be studied. Firstly, virtual damping error is studied.
Transforming Equation (8) into modal coordinates, the EOM of ADM is shown in Equation (31).
Firstly, the error ratio r Am i is derived. When the modal damping method is used, it is shown in Equation (32) where the superscript Am and Lm represents the EOM derive by ADM-M and LMM-M respectively. When s i = 0, r Am i = 0; when s i is small, When s i approaches infinity, r Am i → 1 . The error ratio for different damping ratios and different s i is shown in Table 4. As shown in Table 4 and Equation (33), the errors increase with the increase of frequency, and also increase with the damping ratio of the mode; the error ratio can be 10% at the resonance frequency when ζ = 0.05; when the natural frequency of the mode is smaller than the dominant frequency of the earthquake, the error can also be significant in the higher frequency range. Figure 4 shows examples of r Am for modes with different natural frequencies, the parameter setting is the same as Figure 2.
When i s approaches infinity, → 1 Am i r . The error ratio for different damping ratios and different i s is shown in Table 4. As shown in Table 4 and Equation (33), the errors increase with the increase of frequency, and also increase with the damping ratio of the mode; the error ratio can be 10% at the resonance frequency when ζ = 0.05 ; when the natural frequency of the mode is smaller than the dominant frequency of the earthquake, the error can also be significant in the higher frequency range.  As is shown in Figure 4, the shape of  As is shown in Figure 4, the shape of r Am is not related to ω i . r Am increases monotonously with the increase in frequency and finally approach 1 as the frequency approach infinity. The value of r Am in Figure 4 at the resonance points are marked with squares and the value of them are all 10.0%; the value of r Am when ω = 2ω i are all 19.6%, these values are the same as predicted in Table 4 r Am is caused by neglecting support damping, and is one kind of virtual damping error.
The same as LMM-R, in ADM-R, both the virtual damping error and the damping truncation error exist at the same time.
The error ratio of ADM-R when s i = 1 is shown in Equation (34).
The influence of ζ is small, r Ar i s i =1 also achieves its possible maximum value when ω i = √ ω r 1 ω r 2 ; when s i = 2, r Ar i s i =2 increases linearly with ζ, and increases with the increase of η.
When s i approaches infinity, r Ar i s i →∞ approaches 1. Table 5 shows the value of r Ar i for different parameters. In Table 5, the value of r Ar i s i =2 can be larger than r Ar i s i =1 , indicating that for structures having natural frequencies smaller than the dominant frequency of the excitation, the error caused by the excitation can be larger than the error caused by resonance response. The influence of ν is more significant than ζ at resonance frequency, whereas the influence of ζ is significant when s i = 2. Figure 5 shows the features of r Ar , the parameter setting is the same as Figure 2.   r , indicating that for structures having natural frequencies smaller than the dominant frequency of the excitation, the error caused by the excitation can be larger than the error caused by resonance response. The influence of ν is more significant than ζ at resonance frequency, whereas the influence of ζ is significant when = 2 i s . Figure 5 shows the features of Ar r , the parameter setting is the same as Figure 2. As is shown in Figure 5, the main features of r Ar are the combination of r Am and r Rr . r Ar is the joint result of virtual damping error and damping truncation error. The largest value of r Ar when ω = ω i is slightly larger than that of r Rr ; the values of r Ar when ω = 2ω i are all a bit smaller than r Am but around 8% larger than r Rr .

Damping Formulating Error in MSE Problem
In the previous section, the virtual damping error and the damping truncation error that exist in the FRF of the modal participation factors are studied. In this section, a spectral analysis is used to demonstrate how the effect of MSE affects the accuracy of the structural responses is studied. The spatially variate ground motion is simulated as a random process, this method has been used by many researchers to study the effect of MSE.
The elements in the power spectral density (PSD) matrix of the excitation are usually described as follows.
where S e,ij is the cross power spectral density function between the ground displacement x i and x j . γ ij is the coherence function for x i and x j . γ ij , which is the amplitude of the coherence function, usually named as lagged coherence [16] or incoherence effect [17], is represented by γ l ij (ω) in this study. The phase of γ ij is mainly composed of the wave passage effect, represented as γ w ij (ω), and site response effect describing the phase difference between different soil conditions, represented as γ s ij (ω). Kiureghian [17] proposed the following formula for describing γ ij .
The phases of the coherence function including the wave passage effect and the site-response effect are shown as follows.
where T ij = η ij t ij , describing the time for the ground motion to travel from support i to j; t ij = d ij /v p is the time delay parameter, d ij is the distance between support i and j, v p is the apparent propagation velocity of the earthquake; η ij is a coefficient describing the incident angle of the ground motion, −1 ≤ η ij ≤ 1; T ij is a universal coefficient describing the influence of both v p and η ij ; θ i (ω) is the phase of the local soil filter at support i. Then, the PSD of the absolute displacement responses in modal coordinates is derived. To better illustrate how different factors of MSE affect the structural responses, the soil conditions at different structure supports are assumed to be the same. The PSD of the concerned responses can be derived as follows [1].
is the spectral density function of the ground, S y i is the PSD of the concerned modes y i ; q ij are constant coefficients.
When MSE is not considered, the PSDŜ y i of y i for the corresponding uniform excitation problem is.
Ripples with decreasing amplitude can be observed in Equation (39) because of the decaying parameter γ l jk of the incoherence effect and the cosine function brought by the wave passage effect.
With the help of these ripples, S y i (ω) vibrates around S g f y i 2 n s j=1 q ij 2 . q ij can either be positive or negative, S y i (ω) can be significantly larger or smaller thanŜ y i . When the ridges of the ripples encounter the resonance peak, the amplitude of the resonance response will increase. As a result, the absolute value of the damping truncation error will increase. It is not easy to decide whether the influence of MSE has a positive or negative effect on reducing the virtual damping error and the damping truncation error. The influence of wave passage effect including the influence of the structural span, the incident angle and the apparent propagation velocity of the ground motion is controlled by T jk . When the value of T jk increases, the frequency of the cosine function increases, and the ripples are compressed inside the envelope formed by the incoherence effect. The peaks and troughs of the ripple interact with the resonant peaks of the structure and the energy peaks of ground motion, changing the shape of S y i (ω).
The coherence of the excitation decrease with the increase of the frequency. Thus, the incoherence effect controls the decreasing rate of the ripples' envelope. When the coherence of the excitation approach zero, S y i (ω) approaches S g f y i 2 n s j=1 q ij 2 , in this frequency range the higher modes that are not fully excited in the uniform excitation condition are significantly increased. At the same time, the errors at the resonance peaks would increase.
In the EOM of a structure, the magnitude of the virtual damping error depends on the magnitude of virtual damping in the structural model and the excitation amplitude in absolute coordinates.
Compared with the uniform excitation problem, in the MSE problem, the errors induced by the virtual damping and the damping truncation error have the following new features: 1.
The virtual damping error aroused by less compensated support damping would increase when the amplitude of the excitations in absolute coordinates are larger than their amplitude in relative coordinates.

2.
The interaction of ripples with the resonance peak of the structure can either increase or decrease the damping truncation error and the virtual damping error. 3.
The increase of the energy of higher modes due to the incoherence effect can either increase or decrease the damping truncation error and the virtual damping error. 4.
The requirement of considering more modes in the MSE problem can increase the damping truncation error.

Numerical Simulation
In this section, three cases are used to demonstrate the features of the virtual damping error and the damping truncation error and to show the influence of MSE on these two errors. The EOM of the structures are built with eight modeling methods namely LMM-R/M, RMM-R/M, ADM-R/M, and UN-R/M.

Cascaded Two Supports Multiple Degrees of Freedom Structure
Firstly, cascaded two supports multiple degrees of freedom structures are used to demonstrate the errors exist in the modal frequency response function.
In this part, concentrated multiple degrees of freedom mass-spring-damper systems are used in the numerical simulation to show the features of the virtual damping error and the damping truncation error in the FRF in modal coordinates. The error prediction formulas derived in the previous sections are testified. A sketch drawing of the structure is shown in Figure 6. In this model, n masses and n+1 springs are used. Firstly, cascaded two supports multiple degrees of freedom structures are used to demonstrate the errors exist in the modal frequency response function.
In this part, concentrated multiple degrees of freedom mass-spring-damper systems are used in the numerical simulation to show the features of the virtual damping error and the damping truncation error in the FRF in modal coordinates. The error prediction formulas derived in the previous sections are testified. A sketch drawing of the structure is shown in Figure 6. In this model, n masses and n+1 springs are used. In Figure 6, = 1 m kg, the stiffnesses of the springs are 1 N/m. Rayleigh damping matrix is calculated with the first and the penultimate mode. The damping ratio of the structure is assumed to be 0.05. In the equation of motions derived in the previous section, the excitation coefficients are exactly the same. As a result, when the modal frequency response functions are studied, the form of excitation can only amplify the absolute value of the frequency response function but not the shape of it along the frequency axis. In this simulation, the excitations at the supports are assumed to be identical.
Different numbers of masses are used to study the influence of ν on the damping truncation error. The structure parameters of structures having different numbers of masses are shown in Table   6. For each structure, the mode with the smallest damping ratio is used to calculate μ and ν . These modes are also used for the calculation of the errors in the following context.  In Figure 6, m = 1 kg, the stiffnesses of the springs are 1 N/m. Rayleigh damping matrix is calculated with the first and the penultimate mode. The damping ratio of the structure is assumed to be 0.05. In the equation of motions derived in the previous section, the excitation coefficients are exactly the same. As a result, when the modal frequency response functions are studied, the form of excitation can only amplify the absolute value of the frequency response function but not the shape of it along the frequency axis. In this simulation, the excitations at the supports are assumed to be identical.
Different numbers of masses are used to study the influence of ν on the damping truncation error. The structure parameters of structures having different numbers of masses are shown in Table 6. For each structure, the mode with the smallest damping ratio is used to calculate µ and ν. These modes are also used for the calculation of the errors in the following context. The Fundamental frequency of the structures decreases with the increase of the mass number; whereas µ and ν increases with it. Figure 7 shows the absolute errors and error ratios in the modal FRF of the 4th mode, when n = 20, and when LMM-R and RMM-R are use. The Rayleigh damping matrices are calculated with the 1st and the 19th mode. As shown in Figure 7 a, the input of this frequency response function is acceleration. The value of both the dashed line and the dotted line close to 0 Hz is significant showing the influence of the second term at the right-hand side of Equation (15). The error in this part is the main component of the virtual damping error in LMM-R. One peak of the error and error ratio appears at the resonance point. The shape of the dotted lines around the resonance point in Figure 7 a,b are similar but with small differences, showing that the error in LMM-R near the resonance point contains little virtual damping error and is mainly composed of the damping truncation error. The varying trend of the error ratio is the same as predicted in the previous section. The largest value of the error ratio can be as large as 94% at the resonance point. Figure 8 shows the errors and error ratios of ADM-R and ADM-M when n=20. As shown in Figure 7a, the input of this frequency response function is acceleration. The value of both the dashed line and the dotted line close to 0 Hz is significant showing the influence of the second term at the right-hand side of Equation (15). The error in this part is the main component of the virtual damping error in LMM-R. One peak of the error and error ratio appears at the resonance point. The shape of the dotted lines around the resonance point in Figure 7a,b are similar but with small differences, showing that the error in LMM-R near the resonance point contains little virtual damping error and is mainly composed of the damping truncation error. The varying trend of the error ratio is the same as predicted in the previous section. The largest value of the error ratio can be as large as 94% at the resonance point. Figure 8 shows the errors and error ratios of ADM-R and ADM-M when n = 20.
point. The shape of the dotted lines around the resonance point in Figure 7 a,b are similar but with small differences, showing that the error in LMM-R near the resonance point contains little virtual damping error and is mainly composed of the damping truncation error. The varying trend of the error ratio is the same as predicted in the previous section. The largest value of the error ratio can be as large as 94% at the resonance point. Figure 8 shows the errors and error ratios of ADM-R and ADM-M when n=20. As shown in Figure 8a, the varying trend of the dotted line can be approximated by a linear function which is predicted in Equation (33). The peak value of the dashed line, which is the absolute value of the error, also appears at the resonance point. In Figure 8b, the behavior of the error ratio is more complex than the other three models, multiple peaks appear in the dashed line; but the largest error ratio appears at the resonance point. In the frequency range ω > ω 4 , though the error ratio increases with the increase of frequency, the absolute value of the error decreases with frequency. Comparing Figure 8 with Figure 7, the errors in LMM-R and ADM-R are the joint effect of the virtual damping error and the damping truncation error. Neglecting support damping in ADM-M will induce the virtual damping error. The form of the virtual damping error in ADM-R and the damping truncation error in RMM-R are different. The virtual damping error in ADM-M can influence a wider frequency band and can arise a larger error in the energy of response; whereas the damping truncation error can cause large errors at resonance point.
In Figure 9a,b, the simulated value of the error ratios are compared with the theoretical results introduced in the previous section. Four structures with different numbers of masses are compared. The error ratios are measured at the resonance frequencies of the mode with the smallest damping ratio. The symbolsr Lr ,r Rr ,r Am andr Ar are used, these symbols represent the estimated value of r Lr , r Rr , r Am and r Ar . In Figure 9a, the error ratios in both the LMM-R and RMM-R increase with the increase of n, which reflects the increase in ν. The results show the introduced formula performs well in approximating the error ratio. The error ratios in LMM-R and RMM-R are close to each other but have small differences. ratio. The symbols ˆL r r , ˆR r r , ˆA m r and ˆA r r are used, these symbols represent the estimated value of  Figure 9a, the error ratios in both the LMM-R and RMM-R increase with the increase of n , which reflects the increase in ν . The results show the introduced formula performs well in approximating the error ratio. The error ratios in LMM-R and RMM-R are close to each other but have small differences. In Figure 9b, as stated in the previous section, the error ratio Am r should not have related to the size of the structure. However, the simulation results show Am r changes slightly with the variation of n. For its reason, the error ratio of the numerical simulation is calculated at the resonance peak which is not the natural frequency of the undamped mode but close to it; the FRF changes rapidly in the vicinity of the resonance point contribute to the small variations. The error ratio Ar r also increases with n. The estimated error ratios show high accuracy compared with the simulation results.

4-Span Rigid Frame Bridge
In this part, a 4-span bridge with a total length of 140 m is used to show the influence of MSE on the virtual damping error and the damping truncation error. The data of the bridge can be found in In Figure 9b, as stated in the previous section, the error ratio r Am should not have related to the size of the structure. However, the simulation results show r Am changes slightly with the variation of n. For its reason, the error ratio of the numerical simulation is calculated at the resonance peak which is not the natural frequency of the undamped mode but close to it; the FRF changes rapidly in the vicinity of the resonance point contribute to the small variations. The error ratio r Ar also increases with n. The estimated error ratios show high accuracy compared with the simulation results.

4-Span Rigid Frame Bridge
In this part, a 4-span bridge with a total length of 140 m is used to show the influence of MSE on the virtual damping error and the damping truncation error. The data of the bridge can be found in the reference [2]. Some of the data, including the span and height of the bridge, some of the node number and element number is shown in Figure 10. The frequency response function of the bridge beyond the 8th mode decreases rapidly and the energy of the ground motion is usually small after 10 Hz. Thus, in the analysis, the 1st mode and the 8th mode is used to calculate the Rayleigh damping matrix, and the damping ratio is assumed to be 0.04. A two-dimensional beam element is used as the element of this bridge.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 18 of 27 the reference [2 ]. Some of the data, including the span and height of the bridge, some of the node number and element number is shown in Figure 10. The frequency response function of the bridge beyond the 8th mode decreases rapidly and the energy of the ground motion is usually small after 10 Hz. Thus, in the analysis, the 1st mode and the 8th mode is used to calculate the Rayleigh damping matrix, and the damping ratio is assumed to be 0.04. A two-dimensional beam element is used as the element of this bridge. Spectral analysis is used in this paper to show the statistical feature of the damping formulation problem. The spatially variate ground motions are simulated with stationary random processes.
As is shown in the previous section, LMM-M approaches RMM-M when the added large masses approach infinity. However, in LMM, a limited large number is practically assigned to the large masses instead of infinity. This process induces rigid body modes and modes with very low natural frequencies into the model. These modes do not actually exist in the real structure. If these modes were calculated in the responses, the responses will be overestimated. To study the error in LMM, a trick is indispensable. In this paper, the limit of FRF when the large masses approach infinity is used instead of the practical LMM. And ILMM is used to refer to this ideal LMM. Spectral analysis is used in this paper to show the statistical feature of the damping formulation problem. The spatially variate ground motions are simulated with stationary random processes.
As is shown in the previous section, LMM-M approaches RMM-M when the added large masses approach infinity. However, in LMM, a limited large number is practically assigned to the large masses instead of infinity. This process induces rigid body modes and modes with very low natural frequencies into the model. These modes do not actually exist in the real structure. If these modes were calculated in the responses, the responses will be overestimated. To study the error in LMM, a trick is indispensable. In this paper, the limit of FRF when the large masses approach infinity is used instead of the practical LMM. And ILMM is used to refer to this ideal LMM.
In this study, the soil conditions of the supports are assumed to be the same. The modified Kanai-Tajimi spectrum [1], which is also called the Clough-Penzien spectrum, is used as the spectral density function of the ground motion. The expression of this spectrum is shown in (41).
Firm soil conditions are considered and the parameters are ω g =15 rad, ω f = 1.5 rad, ξ f = 0.6, ξ g = 0.6 Harichandran et al. [18]. S 0 represents the energy intensity of the ground motion. Peak ground acceleration of 0.5 g is used for determining the intensity of the ground motions. The mean peak value of the ground motion and the structural responses are calculated with the method introduced by Der Kiureghian [19]. The duration of the earthquake is assumed to be 20 s.
The model derived by Luco et al. [20] is used to describe the incoherence effect in Equation (37), the expression is shown in Equation (42).
where v s is an estimation for the elastic S-wave velocity in the random medium, λ x = µ √ R/r 0 , R is the distance in the random medium traveled by the wave; r 0 is the scale length of the random inhomogeneities along the path; µ 2 is the measure of the relative variation of the elastic properties in the medium. The value λ x /v s = 2.5 × 10 −4 is recommended in the reference. [15]. The apparent propagation velocity of the ground motion is assumed to be 800 m/s. Existing researches have shown that the apparent propagation velocity is usually larger than 800 m/s. The response PSD derived from different modelling methods are compared. The mean peak responses are used to measure the amplitude of responses.
A new error ratio is defined to measure the error of mean peak responses.
where p stands for the mean peak responses; Method can be replaced by the name of different methods; the method RMM-M is used as the reference method.
The error ratios of all the displacement and force responses when using ILMM are smaller than 10 −8 and are negligible, testifying to the conclusion that LMM approaches RMM when the large masses approaches infinity. Table 7 listed the error ratio of several concerned responses. In Table 7, the responses are represented as D n 6 ,2 and F e 3 ,5 , the former one represents the vertical displacement of node 6; n in the subscript represents the node number; the latter stands for the shear force at the right end of element 3; e in the subscript represents the element number. The key nodes and elements are marked in Figure 10. The numbering of the degree of freedoms in the beam element is also shown in this figure. The numbering of the degree of freedoms in a node is the same as the left node in an element.
As is shown in Table 7, the error ratio of using the Rayleigh damping method in calculating the mean peak response of D n 3 ,2 and F e 3 ,5 are relatively small compared with the other responses. Figure 11 shows the PSD of D n 3 ,2 derived with ILMM and RMM.
The error ratios of all the displacement and force responses when using ILMM are smaller than −8 10 and are negligible, testifying to the conclusion that LMM approaches RMM when the large masses approaches infinity. Table 7 listed the error ratio of several concerned responses. In Table 7, the responses are represented as 6 ,2 n D and 3 ,5 e F , the former one represents the vertical displacement of node 6; n in the subscript represents the node number; the latter stands for the shear force at the right end of element 3; e in the subscript represents the element number. The key nodes and elements are marked in Figure 10. The numbering of the degree of freedoms in the beam element is also shown in this figure. The numbering of the degree of freedoms in a node is the same as the left node in an element. As is shown in Table 7, the error ratio of using the Rayleigh damping method in calculating the mean peak response of 3 ,2 n D and 3 ,5 e F are relatively small compared with the other responses. Figure   11 shows the PSD of 3 ,2 n D derived with ILMM and RMM. As is shown in Figure 11, the responses of ILMM-M and RMM-M overlap with each other testifying that these two modelling methods are equivalent. The response of RMM-R fit very well with RMM-M and starts to lose accuracy near the resonance peak of the 2nd mode, showing the influence of the damping truncation error. This phenomenon is also predicted in the previous derivation. The response of ILMM-R firstly shows an error at around 0.5 Hz, which attributes to the virtual damping error; the response of ILMM-R fit very well with the result of RMM-R showing that ILMM-R has both the virtual damping error and the damping truncation error. The response of D n 3 ,2 is dominated by the responses before the natural frequency of the 1st mode, thus the influence of the damping truncation error is relatively small.
In Table 7, D n 6 ,2 is obviously larger than the other responses. The response PSD of D n 6 ,2 derived with the modelling methods RMM and UN are shown in Figure 12 D is dominated by the responses before the natural frequency of the 1st mode, thus the influence of the damping truncation error is relatively small. In Table 7, 6 ,2 n D is obviously larger than the other responses. The response PSD of 6 ,2 n D derived with the modelling methods RMM and UN are shown in Figure 12. In Figure 12, the response PSD of 6 ,2 n D at the 2nd mode and the 3rd mode is larger relative to the 1st mode, thus the influence of the damping truncation error on 6 ,2 n D is larger than its influence on 3 ,2 n D . Compare the response of RMM and UN, under the influence of MSE the responses of the 2nd mode and the 3rd mode are obviously increased, which largely increases the damping truncation error. In Table 7, the error of ADM-M can either be larger or smaller than the error ratio of RMM-R and ILMM-R showing that the influence of using the Rayleigh damping method can be larger than neglecting the support damping. The error ratio of ADM-R is a joint result of the virtual damping error and the damping truncation error.

Long Span Cable-Stayed Bridge
In the following context, the error in the modelling methods of a cable-stayed bridge is studied. The same analysis as the rigid frame bridge is carried out to this cable-stayed bridge. The bridge is introduced by Dyke et al. [21] and is used as a bench-mark bridge for control methods. The bridge is the Bill Emerson Memorial Bridge spanning the Mississippi River. A simplified drawing is shown in Figure 13. The detailed information of the bridge can be found in Dyke et al. [21]'s article and the data of the bridge can be found on the website. [22]  In Figure 12, the response PSD of D n 6 ,2 at the 2nd mode and the 3rd mode is larger relative to the 1st mode, thus the influence of the damping truncation error on D n 6 ,2 is larger than its influence on D n 3 ,2 . Compare the response of RMM and UN, under the influence of MSE the responses of the 2nd mode and the 3rd mode are obviously increased, which largely increases the damping truncation error.
In Table 7, the error of ADM-M can either be larger or smaller than the error ratio of RMM-R and ILMM-R showing that the influence of using the Rayleigh damping method can be larger than neglecting the support damping. The error ratio of ADM-R is a joint result of the virtual damping error and the damping truncation error.

Long Span Cable-Stayed Bridge
In the following context, the error in the modelling methods of a cable-stayed bridge is studied. The same analysis as the rigid frame bridge is carried out to this cable-stayed bridge. The bridge is introduced by Dyke et al. [21] and is used as a bench-mark bridge for control methods. The bridge is the Bill Emerson Memorial Bridge spanning the Mississippi River. A simplified drawing is shown in Figure 13. The detailed information of the bridge can be found in Dyke et al. [21]'s article and the data of the bridge can be found on the website [22]. The model of the bridge is previously built for uniform excitation. In this paper, the model is rebuilt with the 'hinged' configuration as described in the files of the model. In the 'hinged' configuration, the deck is connected with the towers with shock transmission devices, which can be seen as rigid connections, the data can also be found in the reference. The stiffness matrix of the support is derived from the element stiffness matrices. The damping ratio of the bridge is assumed to be 0.03. The 1st and the 23rd mode are used for the formulation of the Rayleigh damping matrix. The responses under the excitation of ground motions from the longitudinal direction are studied. The soil conditions of the Piers are all assumed to be firm soil, and the energy intensity of the ground motion is decided with the same method as the rigid frame bridge. Table 8 listed the error ratios of the mean peak value of several key responses. The model of the bridge (see Supplementary Materials) is previously built for uniform excitation. In this paper, the model is rebuilt with the 'hinged' configuration as described in the files of the model. In the 'hinged' configuration, the deck is connected with the towers with shock transmission devices, which can be seen as rigid connections, the data can also be found in the reference. The stiffness matrix of the support is derived from the element stiffness matrices. The damping ratio of the bridge is assumed to be 0.03. The 1st and the 23rd mode are used for the formulation of the Rayleigh damping matrix. The responses under the excitation of ground motions from the longitudinal direction are studied. The soil conditions of the Piers are all assumed to be firm soil, and the energy intensity of the ground motion is decided with the same method as the rigid frame bridge. Table 8 listed the error ratios of the mean peak value of several key responses. In Table 8, generally, the error ratio of the cable-stayed bridge is larger than that of the short span rigid frame bridge. The largest damping truncation error of RMM-R can reach 16.15%. The "hinged" model of the cable-stayed bridge assumes that the deck and the two towers of the bridge are perfectly connected, which will definitely increase the responses of the higher modes. Figure 14 shows the response PSD of F d,3 when UN and RMM are used for modelling. As is shown in Figure 14, the higher modes up to 5 Hz has significant energy compared with the first few modes. When using the Rayleigh damping method, the damping ratio of the modes having natural frequencies higher than ω 2 r are larger than ζ 2 r and the resonance responses of these modes are largely suppressed. That is the reason that the error ratio of ,3 d F are large. In Figure 14 As is shown in Figure 14, the higher modes up to 5 Hz has significant energy compared with the first few modes. When using the Rayleigh damping method, the damping ratio of the modes having natural frequencies higher than ω r 2 are larger than ζ r 2 and the resonance responses of these modes are largely suppressed. That is the reason that the error ratio of F d,3 are large. In Figure 14, the PSD of UN-M model has smaller energy at around 4 Hz; however, the PSD of RMM-M has a significant response. This phenomenon is caused by the effect of MSE. In this case, the effect of MSE increases the damping truncation error.
To reduce the error of RMM-R in calculating the responses of F d,3 , using a larger value of ω r 2 can help. However, when ω r 2 increases, the damping ratio of the other responses will definitely increase. The choice of ω r 2 requires more delicate consideration.
In Figure 14, compared with the PSD of UN-M model, due to the influence of MSE, some of the resonance peaks in the PSD of RMM-M are increased and some are decreased. The errors in the UN-R model and the RMM-R model originated from different modes. Therefore, it is hard to decide whether the effect of MSE would increase the damping truncation error or not.
In Table 8, the error ratio of ADM-M and ADM-R have large values when base forces are calculated. The largest error ratio reaches 503%, which is unacceptable showing that the neglect of support damping could cause huge errors. Though the error ratios in ADM of the other responses are acceptable, ADM is not recommended in the applications of the MSE problem.
In this paper, the influence of the apparent propagation velocity of the ground motion is also studied. The influence of the incident angle has a similar effect with the apparent propagation velocity. Figure 15 shows the variation of the mean peak responses of the mid-tower displacement of Pier 2, D m,2 , with respect to the variation of 1/v p . method v is defined as follows to measure the error sensitivity of p v . Table 9 shows the maximum changes of the error ratios when p v changes from 800 m/s to infinity. As is shown in Table 9, the value of

Rr v and
ILr v are mainly around 1%. The largest value is observed from the largest cable force, which is around 3.5%. As is shown in Figure 15, the results of RMM-M and ILMM-M match well with each other. v p mainly influence the vibration frequencies of the ripples caused by the wave passage effect. The error ratios of the other methods either increase or decrease with the variation of 1/v p . v method is defined as follows to measure the error sensitivity of v p . Table 9 shows the maximum changes of the error ratios when v p changes from 800 m/s to infinity. As is shown in Table 9, the value of v Rr and v ILr are mainly around 1%. The largest value is observed from the largest cable force, which is around 3.5%.
The influence of MSE on the error ratio is also studied. Three variables are defined as follows.
In Equation (46), extreme v p (p Rr − p Rm ) represents the extreme value of the damping truncation error of RMM-R when v p varies from 800 m/s to infinity. e Ur is used to measure the error ratio of the damping truncation error when uniform excitation is assumed. r R a is used to measure the relative amplitude of the damping truncation error in the MSE problem and the damping truncation error of the corresponding uniform excitation problem; r R r is used to measure the relative amplitude of error ratio of these two conditions. Figure 16 demonstrates the influence of MSE on the damping truncation error. The error ratios of the key responses are compared in this figure.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 24 of 27 In Equation (46) r is used to measure the relative amplitude of error ratio of these two conditions. Figure 16 demonstrates the influence of MSE on the damping truncation error. The error ratios of the key responses are compared in this figure. the deck-tower connection, the other notation can be found in Table 9. For most responses, the value of R a r is around 1, in other words, the effect of MSE has little influence on the damping truncation error. However, for ,2 t D and ,3 t D , which are the absolute displacements of the tower top, the error in MSE conditions changes to over two times its value in the UN condition; for c F , the error in MSE condition changes more than times six times its value in the UN condition. Looking at the error ratio depicted in Figure 16, it is not hard to conclude that the amplitudes of the error ratios of ,2 R r r shows the changes in error ratio consider or not the MSE, for c F the error ratio in MSE problem is still more than three times that of the uniform excitation problem. For the other responses, the error ratio can increase or decrease by more than 50%. MSE can affect the distribution of error, thus, can significantly increase the damping truncation error.

Conclusions
In this paper, the errors caused by virtual damping error and damping truncation error in the In Figure 16, D t,2 and D t,3 represents the top tower displacement of Pier 2 and Pier 3, respectively; D m,2 and D m,3 are the mid-tower displacements; F dt,2 and F dt,3 are the axial forces of the deck-tower connection, the other notation can be found in Table 9. For most responses, the value of r R a is around 1, in other words, the effect of MSE has little influence on the damping truncation error. However, for D t,2 and D t,3 , which are the absolute displacements of the tower top, the error in MSE conditions changes to over two times its value in the UN condition; for F c , the error in MSE condition changes more than times six times its value in the UN condition. Looking at the error ratio depicted in Figure 16, it is not hard to conclude that the amplitudes of the error ratios of D t,2 and D t,3 have little change except for their sign. The increase in the error is because of the increase in the absolute value of p D t,2 and p D t,3 . r R r shows the changes in error ratio consider or not the MSE, for F c the error ratio in MSE problem is still more than three times that of the uniform excitation problem. For the other responses, the error ratio can increase or decrease by more than 50%. MSE can affect the distribution of error, thus, can significantly increase the damping truncation error.

Conclusions
In this paper, the errors caused by virtual damping error and damping truncation error in the three modeling methods, the large mass method (LMM), the relative motion method (RMM), and the absolute displacement method (ADM), widely used in analyzing multiple support excitation (MSE) problems, are studied.
The term virtual damping is used to define the wrongly induced damping forces generated between the reference coordinates and the structure in the process of damping matrix formulation. The virtual damping error can be eliminated with the formula proposed in this paper. A new Rayleigh damping formula is proposed for LMM to eliminate the virtual damping error. With modal damping and the new Rayleigh damping formula, an integrated equation of motion (EOM) derived from the converged LMM is proved free from the virtual damping error and equivalent to RMM.
The damping truncation error is defined as the error caused by using Rayleigh damping instead of modal damping in the modelling of structures. The error in the structural frequency response function caused by the virtual damping error and the damping truncation error in the aforementioned three modeling methods are studied and the formulas measuring the largest error ratio are derived. The mechanism of how the effect of MSE affects the virtual damping error and the damping truncation error is studied by spectral analysis. A simplified mass-spring model, a short-span rigid-frame bridge and a long-span cable-stayed bridge are used for numerical simulation to testify the conclusions obtained in the theoretical derivation.
The following conclusions are acquired.

1.
RMM is free from the virtual damping error.

2.
When modal damping is used, the LMM converged to RMM. Using Rayleigh damping in LMM consists of both virtual damping error and the damping truncation error. When the newly-introduced Rayleigh damping equation is used in LMM, the LMM also converged to RMM and the virtual damping error in LMM is eliminated. 3.
The virtual damping error exists in all the structural modes of ADM. In ADM, the error ratio of the virtual damping increase with the frequency. 4.
The largest damping truncation errors in LMM, RMM, and ADM appear at the resonance frequency of the mode having the lowest damping ratio. 5.
For a particular structure, the damping truncation errors in LMM, RMM and ADM increase with the parameter ν, which is the quotient of the two frequencies used in formulating the Rayleigh damping matrix. 6.
The effect of MSE changes the energy distribution among modes at the same time changes the distribution of error. Compared with the corresponding uniform excitation problem, the effect of MSE can either increase or decrease the virtual damping error and the damping truncation error in the model. The damping truncation error in the MSE problem can be significantly larger than the error in the corresponding uniform excitation problem.

Conflicts of Interest:
The authors declare no conflict of interest.