An explicit structure-dependent algorithm for pseudodynamic testing

Engineering Structures 46 (2013) 511–525

Contents lists available at SciVerse ScienceDirect

Engineering Structures journal homepage: www.elsevier.com/locate/engstruct

An explicit structure-dependent algorithm for pseudodynamic testing Shuenn-Yih Chang ⇑ Department of Civil Engineering, National Taipei University of Technology, NTUT Box 2653, Taipei 106, Taiwan, ROC

a r t i c l e

i n f o

Article history: Received 6 September 2011 Revised 14 August 2012 Accepted 14 August 2012 Available online 26 September 2012 Keywords: Pseudodynamic testing Structure-dependent Unconditional stability Accuracy Overshooting

a b s t r a c t A novel explicit pseudodynamic algorithm is proposed for the general pseudodynamic testing, where the total response is dominated by low frequency modes while the high frequency responses are of no interest. This is because it can integrate the possibility of unconditional stability and the explicitness of each time step. In fact, it is unconditionally stable for any instantaneous stiffness softening systems, any linear elastic systems, and certain instantaneous stiffness hardening systems that might be experienced in a realistic structure. Hence, stability is not a major concern for this algorithm in practical applications. Meanwhile, it has a second order accuracy and enhanced error propagation properties. As a result, this explicit pseudodynamic algorithm is very promising for the solution of an inertia structural dynamic problem. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction An explicit method [1–7] is generally preferred over an implicit method [8–13] in performing a pseudodynamic test since its implementation is simple when compared to an implicit method. However, there are only a few explicit methods [1–3] in the currently available integration methods and most are implicit methods [1,14–22]. Although many implicit methods are available for the step-by-step integration they are not considered in the early development of the pseudodynamic testing since they involve an iterative procedure and the behavior of a real test specimen is highly path-dependent. Although some implicit pseudodynamic algorithms have been developed for the tests, their implementations become more complex and need some extra hardware [8–13] when compared to explicit pseudodynamic algorithms. In general, two difference equations are needed to form a stepby-step integration method and the coefficients of these two difference equations are specific constants. For this case, it has been shown by Dahlquist [23] that there exists no unconditionally stable explicit method among the multistep methods. This implies that a novel type of integration method is needed so that it can integrate the possibility of unconditional stability and the explicitness of each time step. Learning from the experience in developing the explicit methods with desired numerical dissipation [3,4], an ⇑ Corresponding author. Address: Department of Civil Engineering, National Taipei University of Technology, NTUT Box 2653, No. 1, Section 3, Jungshiau East Road, Taipei 10608, Taiwan, ROC. Tel.: +886 2 2771 2171x2653; fax: +886 2 2781 4518. E-mail address: [email protected] 0141-0296/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.engstruct.2012.08.009

explicit integration method with unconditional stability was developed by Chang on 2002 [24,25], which will be referred as Chang explicit method (CEM). The key step to develop this method is to release the constraint of using constant coefficients for the two difference equations. Hence, the coefficients of the difference equations are no longer constants but allow being functions of initial structural properties and the step size. Later, a similar integration method was proposed by Chen and James on 2008 [26] and is referred as CR explicit method (CRM). Since the coefficients of the difference equations for both CEM and CRM are closely related to the structure analyzed, these two integration methods are structure-dependent. Both methods are shown to be unconditionally stable for linear elastic systems and instantaneous stiffness softening systems while they are conditionally stable for instantaneous stiffness hardening systems. An unusual property is also found for CRM. In fact, it will lead to an overshooting effect for high frequency responses. This implies that the total response might be contaminated or even destroyed by these high frequency responses. In order to improve the stability property for CEM and CRM and to overcome the difficulty from overshooting for CRM, a new explicit method is proposed for the general pseudodynamic testing. Its numerical properties and error propagation properties of the pseudodynamic testing is analytically investigated. Both numerical experiments and actual tests are used to confirm the analytical predictions. 2. Explicit pseudodynamic algorithm In structural dynamics or earthquake engineering, the equation of motion for a single degree of freedom can be expressed as

512

S.-Y. Chang / Engineering Structures 46 (2013) 511–525

€ þ cu_ þ ku ¼ f mu

ð1Þ

where m, c, k and f are the mass, viscous damping coefficient, stiffness and external force, respectively; and u; u_ and ü are the displacement, velocity and acceleration, respectively. An explicit method is proposed for the pseudodynamic testing and it can be expressed as

maiþ1 þ cv iþ1 þ kdiþ1 ¼ fiþ1 diþ1 ¼ di þ ðDtÞv i þ /ðDtÞ2 ai

ð2Þ

v iþ1 ¼ v i þ /ðDtÞai 1

ð3Þ

2 0

1 2

1 þ nX0 þ X

pffiffiffiffiffiffiffiffiffiffiffiffi where X0 = x0(Dt), and x0 ¼ k0 =m is the natural frequency of the system determined from the initial stiffness of k0; and n is a viscous damping ratio. Eq. (2) reveals that the proposed method is explicit since the next step displacement di+1 is directly determined from the second line of this equation, which only involves the data of the (i)th time step. Based on the fundamental theory of structural dynamics, nX0 and X20 can be rewritten in terms of the initial structural properties and the size of integration timeffi step [24,25]. In fact, the relations of pffiffiffiffiffiffiffiffiffiffiffi n = c0/(2mx0) and X0 ¼ k0 =mðDtÞ are substituted into Eq. (3). As a result, the expression of / can be rewritten as

/¼

m mþ

1 ðDtÞc0 2

þ

1 ðDtÞ2 k0 2

ð4Þ

where c0 is introduced to represent an initial viscous damping coefficient and it remains invariant for a whole step-by-step integration procedure. It is found that the determination of the coefficient / using Eq. (4) is much cheaper in computing cost than using Eq. (3) since the expressions of nX0 and X20 will involve solving an eigenvalue problem, which is very time consuming for a multiple degree of freedom system with matrices of larger order. It is worth noting that / is assumed to be invariant in a whole step-by-step integration procedure.

To monitor the stiffness change in the solution of a nonlinear system, the instantaneous degree of nonlinearity di+1 is defined as the ratio of the stiffness at the end of the (i + 1)th time step over the initial stiffness and is

The proposed explicit method (PEM) can be used to perform the step-by-step integration after determining the coefficient /. To formulate the step-by-step integration procedure in a matrix form for the basic analysis, the computing procedure for the (i + 1)th time step is described next. At first, the displacement di+1 can be directly calculated from the second line of Eq. (2) since the acceleration, velocity and displacement at the (i)th time step are available at the beginning of the (i + 1)th time step. Thus, the restoring force ri+1 corresponding to the displacement di+1 can be determined from an assumed force–displacement relation and is ri+1 = ki+1di+1, where ri+1 and ki+1 represent the restoring force and stiffness at the end of the (i + 1)th time step. Next, the velocity at the end of the (i + 1)th time step can be computed by

diþ1  di ¼ Dt

ð5Þ

ð7Þ

This parameter will be used in the subsequent basic analysis and error propagation analysis for PEM. It is clear that di+1 = 1 implies that the instantaneous stiffness at the end of the (i + 1)th time step is equal to the initial stiffness. Whereas, di+1 > 1 represents the case of instantaneous stiffness hardening at the end of the (i + 1)th time step and the case of instantaneous stiffness softening can be represented by 0 < di+1 < 1. It should be mentioned that the instantaneous stiffness hardening or softening adopted herein is used to address that the instantaneous stiffness is larger or less than the initial stiffness. After introducing the instantaneous degree of nonlinearity di+1 and assuming a zero viscous damping ratio, the explicit expression of the amplification matrix Ai+1 is found to be

2

6 Aiþ1 ¼ 4

1

1

/

0

1

/

2 iþ1

X

2 iþ1

X

3

2 iþ1

7 5

ð8Þ

/X

The characteristic equation for this matrix can be derived from the equation of jAi+1  kIj = 0 and is found to be

k3  A1 k2 þ A2 k  A3 ¼ 0

ð9Þ

where k is an eigenvalue of the matrix Ai+1; and the explicit expressions of A1, A2 and A3 are found to be

A1 ¼ 2 

X2iþ1 ; 1 þ 12 X20

A2 ¼ 1;

A3 ¼ 0

ð10Þ

4.1. Stability It is manifested from Eq. (9) that the coefficient of A3 = 0 implies that there is a zero eigenvalue, i.e., k3 = 0. To have a bounded oscillatory response, two stability conditions must be met. One is that there are two complex conjugate eigenvalues, and the other is that jk3j < jk1,2j 6 1. As a result, these two stability conditions lead to

ðA1 Þ2  4A2 6 0;

ð6Þ

0 6 A2 6 1

ð11Þ

After the substitution of Eq. (10) into Eq. (11), the stability conditions for PEM are found to be

0 < X0 6 1

if 2 P diþ1

2 0 < X0 6 Xu0 ¼ pffiffiffiffiffiffiffiffiffiffi if 2 < diþ1

ð12Þ

diþ1 2

ðuÞ

This equation is simply derived from the second and third lines of Eq. (2). Finally, the first line of Eq. (2) is employed to compute the acceleration ai+1. Consequently, this procedure to obtain di+1, vi+1 and ai+1 for the (i + 1)th time step can be written in a recursive matrix form as

Xiþ1 ¼ Aiþ1 Xi þ Liþ1 fiþ1

kiþ1 k0

The numerical properties of PEM can be evaluated from the characteristic equation and will be explored next.

3. Step-by-step computing procedure

v iþ1

4. Numerical properties

diþ1 ¼

where

/¼

 T where Xiþ1 ¼ diþ1 ; ðDtÞv iþ1 ; ðDtÞ2 aiþ1 ; Ai+1 and Li+1 are the amplification matrix and the load vector for the (i + 1)th time step, respectively.

where X0 is used to represent an upper stability limit. Eq. (12) shows that unconditional stability can be achieved for PEM as di+1 6 2 while it can only have conditional stability as di+1 > 2. It is apparent that the conditional stability range highly depends upon the instantaneous degree of nonlinearity. In order to gain an insight into the stability conditions, Eq. (12) is plotted in Fig. 1 for PEM. For comparison, those for CEM and CRM are also shown in Fig. 1. In general, the upper stability limit

513

S.-Y. Chang / Engineering Structures 46 (2013) 511–525

Xiþ1

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4A2 1 ¼ tan  1; A21

niþ1 ¼  lnðA2 Þ 2Xiþ1

ð14Þ

where  niþ1 is a measure of numerical dissipation. The relative period error is often used to measure period distortion and is defined as

Piþ1 ¼

Fig. 1. Variation of upper stability limit with di+1.

decreases from infinity to zero for the three methods except that the starting values to decrease are at 1 for both CEM and CRM and at 2 for PEM. It is also found that the curve for CRM is coincided with that of CEM. This indicates that both methods have the same stability properties. In fact, the unconditional stability range for CEM and CRM are di+1 6 1 and for PEM is di+1 6 2. Meanwhile, the upper stability limit for CEM and CRM are always less than that of PEM in the conditional stability range. Although PEM can only have unconditional stability as di+1 6 2 but not for any instantaneous stiffness hardening systems, this improvement is still very important and is good enough for practical applications. This is because that it is very rare to experience a realistic structure, whose instantaneous stiffness will become twice of the initial stiffness. 4.2. Accuracy To evaluate the period distortion introduced at the (i + 1)th time step, the two complex conjugate eigenvalues of the characteristic equation can be expressed in an exponential form as

k1;2 ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1  A1  j A2  A21 ¼ enXiþ1 jXiþ1 2 4

ð13Þ

pffiffiffiffiffiffiffi  iþ1 ðDtÞ. This expression can be applied where j ¼ 1 and Xiþ1 ¼ x to determine the phase shift of Xiþ1 and the numerical damping ratio of  niþ1 at the (i + 1)th time step. As a result, they can be computed by

T iþ1  T iþ1 xiþ1 ¼ 1  iþ1 T iþ1 x

ð15Þ

 iþ1 and Ti+1 = (2p)/xi+1 are used to represent where T iþ1 ¼ ð2pÞ=x the computed and true periods of the system, respectively. In gen iþ1 are considered as the quantities corresponding to eral,  niþ1 and x ni+1 and xi+1 in a numerical solution. Fig. 2 shows the variations of relative period errors with Dt/T0 as di+1 = 0.5, 1 and 2 for CEM, CRM and PEM. The relative period error generally increases with the increase of Dt/T0 as di+1 is given. In addition, it increases with the decrease of di+1 for a given value of Dt/T0 for each integration method. The curves for CEM and CRM are overlapped together. Hence, they have the same numerical accuracy. In general, PEM shows more period distortion when compared to CEM and CRM for given values of di+1 and Dt/T0. However, this difference in relative period error is insignificant as Dt/ T0 6 0.05. Hence, they have roughly the same accuracy if a time step is chosen to satisfy Dt/T0 6 0.05 for the modes of interest in practical applications. 4.3. Overshoot To show that PEM and CEM do not have the unusual overshooting property found in CRM, a numerical example is examined, where an initial displacement of d0 = 1 with zero initial velocity is considered. The time steps of Dt/T = 10 and 100 are used to obtain the free vibration responses for a linear elastic system. The results are shown in Fig. 3. There is no overshoot in displacement for PEM and CEM, while an overshooting effect is found for CRM and it becomes more significant as Dt/T increases. In fact, it is about 30 times of the initial displacement for Dt/T = 10 and about 300 times for Dt/T = 100. This indicates that a total response may be contaminated or even destroyed by high frequency responses due to the overshooting in displacement. It seems that the linear dependency or linear independency of the eigenvectors as the value of X ? 1 is responsible for the different overshooting phenomena of the three integration methods. In general, each integration method has a zero eigenvalue and two

Fig. 2. Variations of relative period error with Dt/T0.

514

S.-Y. Chang / Engineering Structures 46 (2013) 511–525

Fig. 3. Comparisons of responses obtained from large values of Dt/T for linear elastic system.

complex conjugate eigenvalues for a general value of X. Thus, it has three linearly independent eigenvectors and the convergence of the power of the amplification matrix is guaranteed. In fact, the eigenvector matrices of all the three integration methods are found to be

2

WCEM

1 6 6 ¼6 4 iX

1 iX 2

2 6 ¼4 2 6 ¼4

X

1X 2 4

4

X2

1

1

0

1

1X2 þiX 2 1þ14X2

1X2 iX 2 1þ14X2

X2

X2

1 ffiffiffiffiffiffiffiffiffiffi p

1X2 þiX 2

1þ12X2 1þ12X2

3 2

ð1þ14X2 Þ 7 1 11X2 7 Þ 7 WCRM 2ð 4 5 1þ1X2

1 1þ14X2

3

5. Error propagation properties

7 5 WPEM

1 1X2 iX 2

1 ffiffiffiffiffiffiffiffiffiffi p

1þ12X2 1þ12X2

2

2

X

X

0

3

1 1þ12X2

7 5

ð16Þ

1

It is found that the three eigenvectors of each integration method are different for a general value of X. However, as the value of X ? 1, these three eigenvector matrices reduce to

2

1

1

0

3

1 ¼4 2 X2

1 2 X2

2 3 3 1 1 0 0 i i 0 5 WPEM ¼ 4 1 þ pffiffi2 1  pffiffi2 0 5 1 X2 1 X2

Error propagation analysis of a pseudodynamic algorithm is generally needed [9,29–33]. In the following error propagation analysis, a nonlinear system is considered and the evaluation technique proposed by Chang [34,35] is used. To mathematically model the error propagation procedure in the pseudodynamic testing of a nonlinear system, the following notations are defined. di e

di

a

di

WCEM ¼ 4 iX iX  12 5 WCRM X2 X2 1 2

eigenvalues will become the same and are equal to 1 for CRM as X ? 1. Consequently, the convergence of the power of the amplification matrix is not guaranteed [27,28]. This is the cause to an overshooting effect for CRM. On the other hand, both PEM and CEM can still have three linearly independent eigenvectors as X ? 1 and thus the convergence of power of amplification matrix is guaranteed. As a result, no overshoot was found for PEM and CEM for high frequency modes.

ri

ð17Þ

It is apparent that the two eigenvectors of CRM, which correspond to the two complex conjugate eigenvalues, become the same as X ? 1 and thus CRM no longer has three linearly independent eigenvectors. This is because that two complex conjugate

rei rai edi eri

exact numerically computed displacement at step i without errors. exact displacement at step i, including the effects of errors at previous steps. actual displacement at step i, including the effects of previous errors and errors introduced at the current step. exact numerically computed restoring force at step i without errors. exact restoring force at step i, including the effects of errors at previous steps. actual restoring force at step i, including the effects of previous errors and errors introduced at the current step. displacement error introduced at step i. force error introduced at step i.

515

S.-Y. Chang / Engineering Structures 46 (2013) 511–525

5.1. Error propagation formulation

Using the same evaluation technique proposed in Ref. [34,35], Eq. (25) can be rewritten in the following form

To facilitate the error propagation analysis, the above notations are employed to establish the following relations. a

e

raiþ1

r eiþ1

diþ1 ¼ diþ1 þ ediþ1 ¼

þ

ð18Þ

eriþ1

where the restoring force error eriþ1 can be alternatively expressed as

eriþ1

¼

kiþ1 erd iþ1

ð19Þ

erd iþ1

where is the amount of displacement error corresponding to the restoring force error eriþ1 at the (i + 1)th time step. Using actual displacement and restoring force as shown in Eq. (18), which include the current step errors ediþ1 and eriþ1 , Eq. (6) can be reformulated as

Xeiþ1 ¼ Aiþ1 Xei þ Liþ1 fiþ1 þ Miþ1 edi  Niþ1 erd iþ1 Xeiþ1

where ¼ tors Mi+1 and

2

h

iT ðDtÞ eiþ1 ; ðDtÞ2 aeiþ1 Ni+1 corresponding to edi and e diþ1 ;

3

1 6 7 Miþ1 ¼ 4 0 5; 2 Xiþ1

v

2

0 6 Niþ1 ¼ 4 0

3

2 iþ1

ð20Þ is defined; and the vecerd iþ1 are found to be

7 5

ð21Þ

X

The error propagation equation can be obtained from subtracting Eq. (6) from Eq. (20) and the result is found to be

eiþ1 ¼ Aiþ1 ei þ Miþ1 edi  Niþ1 erdiþ1

ð22Þ

where the error cumulative vector ei+1 for the (i + 1)th time step is defined as eiþ1 ¼ Xeiþ1  Xiþ1 and can be explicitly expressed as

2

3

2 e 3 eð1Þ d  diþ1 iþ1 6 ð2Þ 7 6 iþ1 e 7 7 eiþ1 ¼ 6 4 eiþ1 5 ¼ 4 ðDtÞðv iþ1  v iþ1 Þ 5 ðDtÞ2 ðaeiþ1  aiþ1 Þ eð3Þ iþ1

ð23Þ

ð1Þ

where eiþ1 ¼ ediþ1 is identified and is the cumulative displacement error for the (i + 1)th time step. The matrix Ai+1 and vectors Mi+1 and Ni+1 may vary per time step for a nonlinear system. Hence, the error propagation analysis of a nonlinear system is conducted only for a few consecutive time steps but not for a whole pseudodynamic test procedure. However, it might still give valuable information. For this purpose, the pseudodynamic error, which is propagated from the previous one and two time steps to the current time step, is derived in the following. At first, after the repeated substitutions of ei1 and ei into ei+1 through Eq. (22), the cumulative equation is ðdÞ ðrÞ eiþ1 ¼ eiþ1  eiþ1 ðdÞ eiþ1 ¼ Aiþ1 Ai Mi1 edi2 þ Aiþ1 Mi edi1 þ Miþ1 edi rd rd rd eðrÞ iþ1 ¼ Aiþ1 Ai Ni1 ei1 þ Aiþ1 Ni ei þ Niþ1 eiþ1

ð24Þ

ðdÞ

ðrÞ

ðdÞ

The cumulative error vector ei+1 consists of eiþ1 and eiþ1 , where eiþ1 is the error vector caused by the displacement feedback errors and ðrÞ eiþ1 is that caused by the restoring force feedback errors. After substituting Ak, for k = i and i + 1, and Mk and Nk, for k = i  1, i and i + 1, into Eq. (24), the cumulative displacement error at the (i + 1)th time step is found to be

    rd ediþ1 ¼ Di edi þ Di1 edi1 þ Di2 edi2  Ri erd i þ Ri1 ei1

ð25Þ

where

Di ¼ 1; Di1 ¼

Ri ¼

1 þ 12 X20  X2i 1 þ 12

2 0

X

; Di2 ¼ 1 

2X2i1 þ X2i 1 þ 12

2 0

X

þ

X2i1 X2i 1 þ 12 X20

2

X2i 2X2i1 X2i1 X2i ; R   2 i1 ¼ 2 2 1 þ 12 X0 1 þ 12 X0 1 þ 1 X2 2

0

ð26Þ

ediþ1 ¼ EðdÞ

i1 X

cos½ði  kÞXiþ1 þ aiþ1 edk þ Di2 edi2

k¼i

 EðrÞ

i1 X

sin½ði  kÞXiþ1 þ biþ1 erd k

ð27Þ

k¼i

where ai+1 and bi+1 are phase angles; and the explicit expressions of E(d) and E(r) are found to be vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u  2 u 1 2 u X  X2i 2 iþ1 ðdÞ u   E ¼ t1 þ X2iþ1 1 þ 12 X20  14 X2iþ1 vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   i2 u !2 h 2  u 2 2Xi1 1 þ 12 X20  12 X2i  X2i 1 þ 12 X20  12 X2iþ1 u X ðrÞ i E ¼u þ  2   t 1 þ 12 X20 X2 1 þ 1 X2 1 þ 1 X2  1 X2 iþ1

2

0

2

0

4

ð28Þ

iþ1

where E(d) and E(r) are the error amplification factors for the displacement and restoring force feedback errors from the previous one and two time steps to the cumulative displacement error ediþ1 for the current time step.

5.2. Characteristics of error propagation equation It is manifested from Eq. (27) that the error contribution from the displacement feedback errors edi and edi1 to the cumulative displacement error ediþ1 is dependent upon the error amplification factor of E(d). Whereas, the error contribution from the restoring rd force feedback errors erd i and ei1 to the cumulative displacement error ediþ1 depends upon the error amplification factor of E(r). Thus, the error propagation properties of PEM can be revealed by plotting Eq. (28). The first line of Eq. (28) shows that E(d) is a function of X0, Xi and Xi+1, and hence it varies with di and di+1 only and is independent of di1. On the other hand, the second line of Eq. (28) reveals that E(r) is a function of Xi1 in addition to X0, Xi and Xi+1. Thus, it varies with di1, di and di+1. Although there are many possible combinations of di1, di and di+1, only the cases of di1 = di = di+1 = d = 0.5, 1.0 and 1.5, where each combination is in correspondence to the case of the instantaneous stiffness softening, instantaneous stiffness invariant and instantaneous stiffness hardening, are considered for brevity. Variations of E(d) and E(r) with X0 are shown in Figs. 4 and 5, respectively. In Fig. 4, the error amplification factor E(d) generally increases with the increase of X0 for each curve. This indicates that a high frequency mode will experience more error propagation for each time step when compared to a low frequency mode. It also found that the curve of E(d) for CEM for the case of d = 1.5 tends to infinity as its upper stability limit is approached and it seems to increase linearly for the case of d = 1.0. Whereas, for the case of d = 0.5, it increases gradually from 1 to a larger value and then seemingly levels up as X0 increases from 0 to 10. The characteristic of the curve shown in the case of d = 0.5 for CEM is also found for all the three curves corresponding to the three cases of d = 0.5, 1.0 and 1.5 for PEM. This clearly reveals that PEM has better error propagation properties than for CEM for displacement feedback error. It is also found that the error amplification factor E(d) is increased with the instantaneous degree of nonlinearity for a given X0 for PEM and CEM. Meanwhile, very similar phenomena are also found in Fig. 5 except that the error amplification factor E(r) is no longer starting from 1 but nearly from the origin for both CEM and PEM.

516

S.-Y. Chang / Engineering Structures 46 (2013) 511–525

can be experimentally obtained from the direct stiffness method. In general, this method is to impose a unit displacement for the given degree of freedom by using a servo-hydraulic actuator, and then the restoring forces developed by the specimen in each degree of freedom can be measured by load cells. Hence, these restoring forces give the stiffness coefficients of the correspondent column. This experimental procedure can be repeated until the initial stiffness matrix is achieved. In performing a pseudodynamic test, the second line of Eq. (29) is applied to calculate the displacement vector, and subsequently servo-hydraulic actuators are used to impose the computed displacement vector upon the test specimen. After measuring the restoring force vector through load cells, the velocity vector can be simply obtained from

viþ1 ¼

Fig. 4. Error amplification factor for displacement feedback error.

diþ1  di Dt

ð31Þ

Similar to the manipulation for a single degree of freedom system, this equation is derived from the second and third lines of Eq. (29). Finally, the acceleration vector can be calculated by using the equation of motion, i.e., the first line of Eq. (29). 7. Numerical examples To confirm the numerical properties of PEM for both linear elastic and nonlinear systems some numerical examples are studied. In particular, the enhanced stability property is addressed. In the following numerical study, a two-story shear building is considered and the stiffness for each story is assumed to be in the form of

h i k ¼ k0 1 þ aðDuÞ2

Fig. 5. Error amplification factor for restoring force feedback error.

6. Implementation for MDOF system After obtaining the numerical properties and the error propagation properties of PEM, it is important to further examine the performance of this explicit method in the step-by-step solution of a nonlinear system and in performing a pseudodynamic test. For this purpose, the implementation of PEM for a multiple degree of freedom system is needed since the numerical examples considered and the pseudodynamic tests performed are multiple degree of freedom systems. In general, the formulation of PEM for a multiple degree of freedom system can be expressed as

Maiþ1 þ Cviþ1 þ riþ1 ¼ f iþ1 diþ1 ¼ di þ ðDtÞv i þ UðDtÞ2 ai

ð29Þ

v iþ1 ¼ v i þ UðDtÞai where the restoring force vector may be expressed as ri+1 = Ki+1di+1. The symbol Ki+1 is introduced to represent the stiffness matrix at the end of the (i + 1)th time step. It should be mentioned that the stiffness matrix is not determined during a pseudodynamic test. In correspondence to / for a single degree of freedom system, the coefficient becomes to be a matrix form for a multiple degree of freedom system and is found to be



1 2

U ¼ M þ ðDtÞC0 þ ðDtÞ2 K0

1

M

ð30Þ

where K0 is the initial stiffness matrix; and C0 is the initial damping coefficient matrix and is generally determined from the initial structural properties. Since the coefficient matrix of U must be determined before performing a pseudodynamic test it is of need to obtain the initial stiffness matrix. This initial stiffness matrix

ð32Þ

where k0 is the initial stiffness and Du is a story drift. The structural nonlinearity arises from the story drift if a – 0 is chosen. The case of instantaneous stiffness softening can be simulated if choosing a < 0 while a > 0 can be used to mimic the case of instantaneous stiffness hardening. The lumped mass for the bottom story is taken to be m1 = 103 kg while that for the top story are taken as m2 = 102 kg. In order to confirm the numerical properties of PEM for general nonlinear systems, three structural systems with different stiffness types are considered. The stiffness types include the instantaneous stiffness softening, invariant and hardening. In general, the initial stiffness is 1 2 chosen to be k0 ¼ 108 N=m and k0 ¼ 104 N=m for the bottom and top stories, respectively. The three systems with different stiffness types can be simulated by taking a different coefficient of a for the nonlinear term of each story. System-A a1 = 0

a2 = 0

System-B a1 =  1.15  107 a2 =  1.35  101

System-C a1 = 1.15  107

a2 = 1.35  101

instantaneous stiffness invariant system instantaneous stiffness softening system instantaneous stiffness hardening system

where a1 and a2 are the constant coefficients of the nonlinear stiffness terms for the bottom and top stories, respectively. The natural frequencies of the structure are found to be 10.00 and 316.24 rad/s based on the initial stiffness matrix. This structure

S.-Y. Chang / Engineering Structures 46 (2013) 511–525

517

is excited by a ground acceleration of 10 sin(2t) at the base. For each system, the numerical results obtained from the Newmark explicit method (NEM) with Dt = 0.0001 s are considered as ‘‘exact’’ solutions for comparison.

analytical predictions of relative period errors are consistent with the numerical results shown in Fig. 6. Unconditional stability is indicated by these results for both methods since no numerical ð2Þ explosions occur although the value of x0 ðDtÞ is as large as 9.49.

7.1. System-A

7.2. System-B

Numerical results for System-A are plotted in Fig. 6. It is anticipated that both PEM and CRM will lead to acceptable solutions if the time step is selected based on accuracy consideration only since they are unconditionally stable for a linear elastic system. It is manifested from Fig. 2 that a time step will lead to an accurate integration of the mode of interest if Dt/T0 6 0.05 is met. As a result, the time step of Dt = 0.03 s is chosen to perform the stepby-step integration. This time step leads to a relative period error of 0.75% for CRM and 1.86% for PEM for the first mode. These

The displacement responses for System-B are shown in Fig. 7. Meanwhile, the time histories of relative period error, instantaneous of degree of nonlinearity and upper stability limit for the two modes are shown in Fig. 8. It is of no need to consider stability problem for PEM and CRM since System-B is an instantaneous stiffness softening system. The time step of Dt = 0.03 s is also chosen for the step-by-step integration. It is manifested from Fig. 8a that the first mode can be reliably integrated and Fig. 8b that a significant period distortion is found in the second mode.

Fig. 6. Displacement responses of System-A.

Fig. 7. Displacement responses of System-B.

518

S.-Y. Chang / Engineering Structures 46 (2013) 511–525

Fig. 8. Response time histories of System-B.

Fig. 9. Displacement responses of System-C.

Since reliable solutions are obtained from PEM and CRM with Dt = 0.03 s the total response is dominated by the first mode and the second mode response is negligible. Fig. 8c and d shows that di varies only between 1 and 0.7 for both modes and thus System-B is an instantaneous stiffness softening system. This implies that both CRM and PEM can have unconditional stability in the step-by-step solution of this system. This is confirmed by the numerical results shown in Fig. 7 for both PEM and CRM.

7.3. System-C An instantaneous stiffness hardening system is examined in this example. In order to show the superiority of PEM over CRM in the step-by-step solution of an instantaneous stiffness hardening system, a time step of Dt = 0.03 s is also chosen to carry out the step-by-step integration based on accuracy consideration. Numerical results are plotted in Fig. 9. Although this time step is satisfied

S.-Y. Chang / Engineering Structures 46 (2013) 511–525

519

Fig. 10. Response time histories of System-C.

Fig. 11. Seismic response of bottom story for 2-story structure with simulation errors.

the rough guideline of Dt/T0 6 0.05 to obtain a reliable solution for PEM and CRM, only PEM leads to an acceptable solution while that obtained from CRM blows up very early. This may be explained by

Fig. 10, where the time histories of relative period error, instantaneous degree of nonlinearity and upper stability limit for both ð1Þ ð2Þ modes are plotted. Since Fig. 10c and d shows that di and di

520

S.-Y. Chang / Engineering Structures 46 (2013) 511–525

Fig. 12. Seismic response of top story for 2-story structure with simulation errors.

Fig. 13. Frequency content of pseudodynamic response of each degree of freedom.

are always greater than 1 and less than about 1.2 for both modes, the structure is an instantaneous stiffness hardening system. Hence, CRM can only have conditional stability and the instability

is due to the violation of the upper stability limit. On the other hand, PEM still gives a very reliable solution since it is unconditionally stable for the range of 0 < di 6 2. Fig. 10a implies that the

521

S.-Y. Chang / Engineering Structures 46 (2013) 511–525

dominant first mode can be reliably integrated since the relative period error is only about 1.85% for PEM, and thus an acceptable solution is obtained. 8. Computer simulations of pseudodynamic testing



1 501

;

 /2 ¼

2003 1

Mode

1st mode

2nd mode

x (rad/s) X = x(Dt)

5.00 0.100 1.001 1.001 0.100 0.100

223.83 4.477 2.453 1.354 4.480 1.823

(d)

E

To numerically illustrate that PEM has much better error propagation properties than those of CEM, computer simulations of these two pseudodynamic algorithms are conducted next. In these simulations, the displacement error introduced into the imposed displacement in each time step is considered as a random variable, with a distribution assumed to be a truncated normal distribution. Simulation details can be found in Ref. [24] and will not be elaborated here. In order to simulate a properly adjusted test, the mean value of the truncated normal distribution is taken to be zero. In addition, its standard deviation is assumed to be 1/3 of the tolerance limit, which is taken to be 0.03 mm in each numerical simulation. This implies that the random error simulated in each time step might be as large as 0.03 mm, which is about the maximum step error that may actually occur in performing the pseudodynamic test. The structural system used in the previous section is adopted here with some modifications for computer simulations. The lumped mass for the bottom and top stories are taken as m1 = 1 kg and m2 = 4 kg, while the initial stiffness is selected to 1 2 be k0 ¼ 5  104 N=m for the bottom story and k0 ¼ 102 N=m for the top story. The system is assumed to be linear elastic, i.e., a1 = a2 = 0. The natural frequencies of the structural system are found to be 5.00 and 223.83 rad/s and their correspondent mode shapes are

/1 ¼

Table 1 Comparisons of error amplification factors.

ð31Þ

The ground accelerations recorded from 1971 San Fernando earthquake are chosen for the seismic inputs, where the peak ground acceleration is scaled to 0.5 g. Both PEM and CEM are used to compute the responses for the cases with and without simulation errors. The time step of Dt = 0.02 s is used in all the computations and the results are plotted in Figs. 11 and 12 for the bottom and top stories, respectively. Each figure contains three plots, which show the time histories of the displacement responses, simulated step displacement errors and the cumulative displacement errors from the top plot to the bottom plot. Furthermore, the frequency content of the displacement response of each degree of freedom for the simulation results obtained from both pseudodynamic algorithms are also computed and shown in Fig. 13. In addition, error amplification factors for both modes of these two methods are listed in Table 1 for the subsequent discussions and comparisons. In Fig. 11, the bottom story response obtained from CEM is highly plagued by simulation errors while that obtained from PEM is only slightly contaminated by simulation errors. On the other hand, Fig. 12 reveals that the top story responses for both methods are almost unaffected by the simulation errors since the two curves with and without simulation errors are almost overlapped. These different phenomena can be explained by two considerations. First, the top story response is dominated by the first mode only. Whereas, in addition to the contribution from the first mode, the contribution from the second mode to the bottom story response is considerable. This difference in response contribution for the top and bottom stories is because that the response contribution from the second mode to the bottom story is much larger than that to the top story. In fact, it is in a ratio of – 2003, which is manifested from the second mode shape /2 in Eq. (31). Second, error amplification factors, both E(d) and E(r), were small and roughly equal for the first mode of the two pseudodynamic algo-

(r)

E

CEM PEM CEM PEM

rithms while they became quite different for the second mode as shown in Table 1. Therefore, very accurate top story responses can be obtained for both algorithms since the first mode is integrated accurately, the response contribution from the second mode is negligible, and error amplification factors are very small for the first mode of the two algorithms. Whereas, the quite differences in error amplification factors of the second mode between PEM and CEM are responsible for the different bottom story responses. This is because that the error amplification factors for the second mode of CEM are E(d) = 2.453 and E(r) = 4.480, which are larger than E(d) = 1.354 and E(r) = 1.823 for PEM. Thus, the bottom story response obtained from CEM is less accurate than that of PEM due to errors associated with the second mode response. Other phenomena are also found in Figs. 11–13. In Fig. 11, it shows that the simulated step errors introduced for each time step are almost the same for the two algorithms at both the top and bottom stories. The cumulative errors for CEM at the bottom story are much larger than those for PEM while they are likely the same for the top story. This is consistent with the displacement responses shown in Figs. 11a and 12a. It is manifested from Fig. 13a that the first mode response of the bottom story obtained from CEM is almost not contaminated by simulation errors since the numerical results with and without simulation errors have roughly the same response amplitude of 0.31 mm for the first mode. Whereas, the second mode response of the bottom story obtained from CEM is badly contaminated by simulation errors since the response amplitude of 1.09 mm is found for the second mode for the numerical results with simulation errors, which is much larger than that of 0.41 mm for those without simulation errors. Hence, the total response is almost destroyed. Very similar phenomena are also found in Fig. 13c for the bottom story responses obtained from PEM, except that the error propagation effect for the second mode response is less significant than for CEM. In fact, the response amplitude of 0.23 mm is found for the second mode for the numerical results with simulation errors and is only slightly larger than that of 0.13 mm for those without simulation errors. As a result, the total response seems acceptable. On the other hand, Fig. 13b and d shows that the displacement responses obtained from the two algorithms at the top story almost originated from the first mode only. They also reveal that the first mode response is almost unaffected by the simulation errors since the two curves are almost coincided. This figure also shows that the period of the first mode was almost not distorted for the two algorithms because the value of Dt/T is only 0.032 for this mode. However, the period of the second mode with Dt/T = 0.45 is elongated from 37.31 Hz to 18.31 Hz and 11.78 Hz for CEM and PEM, respectively. 9. Pseudodynamic tests In order to verify the feasibility of PEM and its superior properties over CRM, a series of pseudodynamic tests were conducted. In fact, a 2-degree-of-freedom system was designed and fabricated. Several steel beams with the cross section of H100  100  6  8 and a length of 2.2 m were adopted for the tests. The test setup is shown in Fig. 14 where the cantilever beam is loaded by 2

522

S.-Y. Chang / Engineering Structures 46 (2013) 511–525

actuators in parallel for simulating a 2-degrees-of-freedom system. To avoid the problem of local buckling or instability, the steel beam was loaded in its minor axis. Before performing an actual pseudodynamic test, the initial stiffness matrix of the specimen can be experimentally measured and was found to be

" K0 ¼

5

6

8:3  10

1:7  10

1:7  106

4:4  106

# ð32Þ

where the unit for each element is in kN/mm. The measured offdiagonal term kij is not exactly equal to kji for i – j. However, an average value is used in order to have a symmetric stiffness matrix, which will be used to compute coefficient matrices. On the other hand, two lumped mass matrices were specified as follows so that two different systems S1 and S2 can be simulated.

" M¼ " M¼

5  104 0 5  104 0

0

#

for S1 102 # 0 for S2 1

ð33Þ

where the unit for each element is in kg. As a result, the initial natural frequencies of S1 are found to be 1.86 and 209.79 rad/s for the first and second modes, respectively, while those for S2 are 1.86 and 2097.62 rad/s. It is clear that both S1 and S2 are designed to have a relatively high second mode. The specimen is subjected to 1971 San Fernando earthquake with a peak ground acceleration (PGA) of 0.4 g for linear elastic tests while for inelastic tests PGA is scaled to 1.0 g. Both CRM and PEM with the time step of Dt = 0.04 s are employed to carry out the step-by-step integration. Pseudodynamic test results for S1 are shown in Fig. 15 for linear elastic tests and in Fig. 16 for

Fig. 14. Test setup for a 2-degree of freedom system.

Fig. 15. Pseudodynamic responses to San Fernando earthquake for S1 (PGA = 0.4 g).

S.-Y. Chang / Engineering Structures 46 (2013) 511–525

523

Fig. 16. Pseudodynamic responses to San Fernando earthquake for S1 (PGA = 1.0 g).

Fig. 17. Pseudodynamic responses to San Fernando earthquake for S2 (PGA = 0.4 g).

nonlinear tests, while those for S2 are shown in Fig. 17 for linear tests and in Fig. 18 for nonlinear tests. Pseudodynamic responses obtained from CEM with a time step of Dt = 0.02 s are considered as ‘‘correct’’ solutions for comparisons. For this time step, the values of D t/T0 for the first and second modes are found to be 0.0059 and 0.66, respectively, for S1, while for S2 they are found to be 0.0059 and 6.67 correspondingly. Apparently, the first mode can be very accurately integrated for both linear elastic tests and inelastic tests for S1 and S2, while very significant period distortion is found for the second mode. Since the response contribution from the second mode response to the total response is insignificant for this problem, very reliable results are achieved and thus can be considered as reference solutions. In general, both CRM and PEM with Dt = 0.04 s can provide an acceptable result for S1 for both linear elastic and inelastic tests as shown in Figs. 15 and 16 although the results obtained from CRM show a little bit more fluctuation in the second mode

response than those obtained from PEM in the top story response for both linear elastic and inelastic tests. On the other hand, for S2, the pseudodynamic results obtained from CRM with Dt = 0.04 s are destroyed by the spurious participation of second mode for both linear elastic and inelastic tests while it seems that reliable pseudodynamic results can still be obtained from PEM with Dt = 0.04 s for both linear elastic and inelastic tests. All these results are manifested from Figs. 17 and 18. These results can be explained next. In general, the first mode response can be very accurately integrated for using Dt = 0.04 s for both CRM and PEM while the second mode response is significantly distorted. Since the contribution from the second mode to the total response is small and thus it is negligible, PEM gives reliable results for both linear elastic and inelastic tests for both S1 and S2. However, an overshooting effect will be found in the second mode for using CRM since the values of Dt/T are as large as 1.34 and 13.35 for the second mode of S1 and S2, respectively. Although the response

524

S.-Y. Chang / Engineering Structures 46 (2013) 511–525

Fig. 18. Pseudodynamic responses to San Fernando earthquake for S2 (PGA = 1.0 g).

contribution from the second mode to the total response is small, the second mode response was significantly amplified due to overshooting for using CRM and thus the total response was contaminated for S1 and was destroyed for S2. Since the value of Dt/T is 1.34 for S1 and 13.35 for S2 in the second mode, the overshooting effect for S2 is more significant than that of S1 in the second mode for using CRM and thus the top story responses obtained from CRM for S1 were contaminated by experimental errors while those for S2 were destroyed by experimental errors. Unconditional stability for linear elastic systems and any instantaneous stiffness softening systems for CEM, CRM and PEM were also indicated by these tests ð2Þ since they all remained stable as the value of X0 was as large as 83.90 for the second mode for both linear elastic and inelastic tests.

10. Conclusions A pseudodynamic algorithm is proposed herein. It can be implemented as simple as an explicit pseudodynamic algorithm since it involves no nonlinear iterations. In addition, the most important enhanced property of this explicit algorithm is that it can have unconditional stability for a certain instantaneous hardening systems in addition to any instantaneous stiffness softening systems and linear elastic systems. In fact, it is unconditionally stable for an instantaneous stiffness hardening system with the instantaneous degree of nonlinearity less than or equal to 2. This improvement allows it to conduct a general pseudodynamic test without considering the stability problem since it is very rare to experience an actual structural system, whose instantaneous degree of nonlinearity is greater than 2. Hence, this pseudodynamic algorithm is much better than Chang explicit method and CR explicit method since they can only have conditional stability for an instantaneous stiffness hardening system. In addition to stability, this algorithm can have comparable accuracy as that of a second-order accurate method, such as the Newmark explicit method and Chang explicit method. It is also found that both the proposed algorithm and Chang explicit method have no overshooting in high frequency modes and thus their performance will be much better than CR explicit method for a system with high frequency modes. Numerical examples are used to confirm the stability and accuracy properties of the proposed explicit algorithm. Furthermore, computer simulations of pseudodynamic tests are used to confirm

the superior error propagation properties of the proposed explicit algorithm in contrast to Chang explicit method. Finally, the actual pseudodynamic tests were performed to confirm the feasibility of the proposed explicit algorithm and its superiority over CR explicit method for a system with high frequency mode. Consequently, this pseudodynamic algorithm is very promising for the solution of a general structural dynamic problem, where the total response is dominated by low frequency modes while the high frequency responses are of no interest. Acknowledgement The author is grateful to acknowledge that this study is financially supported by the National Science Council, Taiwan, ROC, under Grant No. NSC-99-2221-E-027-029. References [1] Newmark NM. A method of computation for structural dynamics. J Eng Mech Div, ASCE 1959;85:67–94. [2] Shing PB, Mahin SA. Elimination of spurious higher-mode response in pseudodynamic tests. Earthq Eng Struct Dynam 1987;15:425–45. [3] Chang SY. Improved numerical dissipation for explicit methods in pseudodynamic tests. Earthq Eng Struct Dynam 1997;26:917–29. [4] Chang SY. The c-function pseudodynamic algorithm. J Earthq Eng 2000;4(3):303–20. [5] Chang SY. Application of the momentum equations of motion to pseudodynamic testing. Philos Trans Roy Soc, Ser A 2001;359(1786):1801–27. [6] Chang SY. An improved on-line dynamic testing method. Eng Struct 2002;24(5):587–96. [7] Chang SY, Tsai KC, Chen KC. Improved time integration for pseudodynamic tests. Earthq Eng Struct Dynam 1998;27:711–30. [8] Nakashima M, Kaminosomo T, Ishida M. Integration techniques for substructure pseudodynamic test. In: Proceeding of fourth US national conference on earthquake engineering, vol. 2; 1990. p. 515–24. [9] Shing PB, Vannan MT, Carter E. Implicit time integration for pseudodynamic tests. Earthq Eng Struct Dynam 1991;20:551–76. [10] Chang SY, Mahin SA. Two new Implicit algorithms of pseudodynamic test methods. J Chin Inst Eng 1993;16(5):651–64. [11] Thewalt CR, Mahin SA. An unconditionally stable hybrid pseudodynamic algorithm. Earthq Eng Struct Dynam 1995;24:723–31. [12] Chung W, Campbell SD. Implicit pseudodynamic algorithm with an event-toevent solution scheme. Eng Struct 2007;29(4):640–8. [13] Hung CC, Yen WM, Lu WT. An unconditionally stable algorithm with the capability of restraining the influence of actuator control errors in hybrid simulation. Eng Struct 2012;42:168–78. [14] Houbolt JC. A recurrence matrix solution for the dynamic response of elastic aircraft. J Aeronaut Sci 1950;17:540–50.

S.-Y. Chang / Engineering Structures 46 (2013) 511–525 [15] Wilson EL, Farhoomand I, Bathe KJ. Nonlinear dynamic analysis of complex structures. Earthq Eng Struct Dynam 1973;1:241–52. [16] Park KC. An improved stiffly stable method for direct integration of nonlinear structural dynamic equations. J Appl Mech 1975;42:464–70. [17] Hilber HM, Hughes TJR, Taylor RL. Improved numerical dissipation for time integration algorithms in structural dynamics. Earthq Eng Struct Dynam 1977;5:283–92. [18] Wood WL, Bossak M, Zienkiewicz OC. An alpha modification of Newmark’s method. Int J Numer Meth Eng 1981;15:1562–6. [19] Chung J, Hulbert GM. A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-a method. J Appl Mech 1993;60(6):371–5. [20] Chang SY. A series of energy conserving algorithms for structural dynamics. J Chin Inst Eng 1996;19(2):219–30. [21] Tamma KK, Zhou X, Sha D. A theory of development and design of generalized integration operators for computational structural dynamics. Int J Numer Meth Eng 2001;50:1619–64. [22] Zhou X, Tamma KK. Design, analysis and synthesis of generalized single step solve and optimal algorithms for structural dynamics. Int J Numer Meth Eng 2004;59:597–668. [23] Dahlquist G. A special stability problem for linear multistep methods. BIT 1963;3:27–43. [24] Chang SY. Explicit pseudodynamic algorithm with unconditional stability. J Eng Mech, ASCE 2002;128(9):935–47.

525

[25] Chang SY. Improved explicit method for structural dynamics. J Eng Mech, ASCE 2007;133(7):748–60. [26] Chen C, Ricles JM, Marullo TM, Mercan O. Real-time hybrid testing using the unconditionally stable explicit CR integration algorithm. Earthq Eng Struct Dynam 2009;38:23–44. [27] Hughes TJR. The finite element method. Englewood Cliffs (NJ, USA): PrenticeHall, Inc; 1987. [28] Belytschko T, Hughes TJR. Computational methods for transient analysis. North-Holland: Elsevier Science Publishers BV; 1983. [29] Peek R, Yi WH. Error analysis for pseudodynamic test method. I: analysis. J Eng Mech, ASCE 1990;116:1618–37. [30] Peek R, Yi WH. Error analysis for pseudodynamic test method. II: application. J Eng Mech, ASCE 1990;116:1638–58. [31] Shing PB, Mahin SA. Cumulative experimental errors in pseudo dynamic tests. Earthq Eng Struct Dynam 1987;15:409–24. [32] Shing PB, Mahin SA. Experimental error effects in pseudodynamic testing. J Eng Mech, ASCE 1990;116:805–21. [33] Shing PB, Manivannan T. On the accuracy of an implicit algorithm for pseudodynamic tests. Earthq Eng Struct Dynam 1990;19:631–51. [34] Chang SY. Nonlinear error propagation analysis for explicit pseudodynamic algorithm. J Eng Mech, ASCE 2003;129(8):841–50. [35] Chang SY. Error propagation in implicit pseudodynamic testing of nonlinear systems. J Eng Mech, ASCE 2005;131(12):1257–69.