# Theory and Implementation of a Two-Step Unconditionally Stable Explicit Integration Algorithm for Vibration Analysis of Structures

Authors: Li Changqing; Jie Junping; Jiang Lizhong; T. Y. Yang
Journal: Shock and Vibration
 (2016)

---

## Abstract

Time history analysis is becoming the routine process to quantify the response of the structure under dynamic loads. In this paper, a novel two-step unconditionally stable explicit integration algorithm, named Unconditional Stable Two-Step Explicit Displacement Method (USTEDM), is proposed for vibration analysis of structure. USTEDM is unconditionally stable, requires low memory, produces no overshoot, and is third order accurate. The accuracy and efficiency of USTEDM are presented and compared with other commonly used integration algorithms. The result shows that the proposed algorithm has superior performance and can be used efficiently in solving vibration response of civil engineering structure.


---

## Body

## 1. Introduction

With the advance of computer hardware, more advanced computer software is being used to solve the response of the structure under dynamic loads. Hence, more advanced numerical methods have been developed to efficiently and accurately predict the dynamic response of the structure. In the past, multiple implicit [1–7] and explicit [8–12] integration algorithms have been developed, whose relationship [13] and design methods are comprehensively analyzed [14–16]. The explicit method has high computational efficiency, while the implicit method has higher stable. The commonly referenced explicit algorithm is Central Difference Method (CDM) (ref), CDM is easy to use, but it is only conditionally stable. Many new explicit methods which are unconditionally stable [17–19] have been developed. Though these methods are unconditionally stable, they require multiple steps which are not as efficient as the CDM. The paper presents a novel 2-step explicit integration method, named Unconditional Stable Two-Step Explicit Displacement Method (USTEDM), that is used to solve the dynamic response of a structure. USTEDM is efficient to use and unconditionally stable. The performance of USTEDM is compared with other explicit algorithms. The result shows that USTEDM has superior performance and can be used efficiently in solving vibration response of civil engineering structures.


## 2. Derivation of the Two-Step Integration Algorithm

Equation (1) shows the generalized dynamic equations of motion of a structure:(1)mat+cvt+kut=ft,where a(t),v(t),u(t) are the acceleration, velocity, and displacement of the structure, respectively. m,c,k are the mass, damping, and stiffness of the structure, respectively. f(t) is external applied force.In this paper, a two-step integration algorithm, as presented in (2), is proposed to solve the dynamic equation as presented in (1):(2)ut+2Δt=fut+Δt,ut,ft+2Δt,ft+Δt,ft,m,c,k,Δt,where ut+2Δt is the displacement at time step t+2Δt and Δt is the integration time step. Similarly, ut+Δt and ut are the displacement at time t+Δt and t, while ft+2Δt,ft+Δt and ft are the force at t+2Δt,t+Δt, and t, respectively.It should be noted that (2) as presented above has physical units, where the displacement could be in the unit of “meter” or “feet,” while the mass can be in the unit of “kg” or “slugs.” To make the equation dimensionless, (2) can be rewritten in form of dimensionless quantity as(3)ut+2Δtut+Δt=ut+2Δtut+Δt1,utut+Δt,ft+2ΔtΔt2mut+Δt,ft+ΔtΔt2mut+Δt,ftΔt2mut+Δt,1,cΔtm,kΔt2m,1.To construct functionut+2Δt/ut+Δt at the right hand of (3), the linear combination of ut+2Δt/ut+Δt with the arguments including 1, (ut/ut+Δt)(ft+2ΔtΔt2/mut+Δt),ft+ΔtΔt2/mut+Δt and ftΔt2/mut+Δt, can be used:(4)r1ut+2Δtut+Δt=r2+r3utut+Δt+r4ft+2ΔtΔt2mut+Δt+r5ft+ΔtΔt2mut+Δt+r6ftΔt2mut+Δt,where rii=1,…,6 are undetermined coefficients. The effects of the dimensionless arguments cΔt/m,kΔt2/m that represent the structure’s material properties to the displacements can be expressed by replacing the coefficients ri(i=1,2,3) with linear combination of constant term 1,cΔt/m,kΔt2/m:(5)r1=l1+l2cΔtm+l3kΔt2m,r2=l4+l5cΔtm+l6kΔt2m,r3=l7+l8cΔtm+l9kΔt2m,where li(i=1,…,9) are coefficients.Substituting (5) into (4) yields(6)l1+l2cΔtm+l3kΔt2mut+2Δtut+Δt=l4+l5cΔtm+l6kΔt2m+l7+l8cΔtm+l9kΔt2mutut+Δt+l10ft+2ΔtΔt2mut+Δt+l11ft+ΔtΔt2mut+Δt+l12ftΔt2mut+Δt,where li(i=1,…,12) are coefficients to be determined and l10=r4,l11=r5,l12=r6.Equation (6) can be simplified to be(7)l1m+l2cΔt+l3kΔt2ut+2Δt=l4m+l5cΔt+l6kΔt2ut+Δt+l7m+l8cΔt+l9kΔt2ut+l10ft+2ΔtΔt2+l11ft+ΔtΔt2+l12ftΔt2.In order to decrease the number of undetermined coefficientsli, (7) can be simplified by dividing l1 at both sides:(8)m+b1cΔt+b2kΔt2ut+2Δt=b3m+b4cΔt+b5kΔt2ut+Δt+b6m+b7cΔt+b8kΔt2ut+b9ft+2ΔtΔt2+b10ft+ΔtΔt2+b11ftΔt2,where bi=(li+1/l1)(i=1,…,11) are undetermined coefficients.Equation (8) is the preliminary form of USTEDM.


## 3. Determination of Coefficients

### 3.1. Accuracy Analysis

The accuracy of USTEDM can be examined using convergence test. Equation (9) shows the Taylor series expansion of the displacement, velocity, and acceleration of the structure at time step of t+2Δt, t+Δt, respectively.(9)ut+Δt=ut+vtdt+12atdt2+16a˙tdt3+⋯,vt+Δt=vt+atdt+12a˙tdt2+16a¨tdt3+⋯,at+Δt=at+a˙tdt+12a¨tdt2+16a⃛tdt3+⋯,ut+2Δt=ut+2vtdt+2atdt2+43a˙tdt3+⋯,vt+2Δt=vt+2atdt+2a˙tdt2+43a¨tdt3+⋯,at+2Δt=at+2a˙tdt+2a¨tdt2+43a⃛tdt3+⋯.Similarly, the dynamic equilibriums att+2Δt,t+Δt,t are expressed in(10)ft+2Δt=mat+2Δt+cvt+2Δt+kut+2Δt,ft+Δt=mat+Δt+cvt+Δt+kut+Δt,ft=mat+cvt+kut.Substituting (9) into (10) yields(11)ft+2Δt=mat+2a˙tdt+⋯+cut+2vtdt+⋯+kut+2vtdt+⋯,ft+Δt=mat+a˙tdt+⋯+cvt+atdt+⋯+kut+vtdt+⋯,ft=mat+cvt+kut.Substituting (9)–(11) into (8) yields(12)m+b1cΔt+b2kΔt2εut+2Δt=m+b1cΔt+b2kΔt2ut+2Δt-b3m+b4cΔt+b5kΔt2ut+Δt-b6m+b7cΔt+b8kΔt2ut-b9ft+2ΔtΔt2-b10ft+ΔtΔt2-b11ftΔt2=d0+d1Δt+d2Δt2+d3Δt3+⋯,where(13)d0=1-b3-b6utm,d1=2-b3vtm+b1-b4-b7utc,d2=2-0.5b3-b9-b10-b11atm+2b1-b4-b9-b10-b11vtc+b2-b5-b8-b9-b10-b11utk,d3=43-16b3-2b9-b10a˙tm+2b1-0.5b4-2b9-b10atc+2b2-b5-2b9-b10vtk,d4=23-124b3-2b9-12b10a¨tm+43b1-16b4-2b9-12b10a˙tc+2b2-12b5-2b9-12b10atk,d5=415-1120b3-43b9-16b10a⃛tm+23b1-124b4-43b9-16b10a¨tc+43b2-16b5-43b9-16b10a˙tk⋯.USTEDM having first-order accuracy means(14)limΔt→0εut+2ΔtΔt=0.Combining (12) and (14) leads to(15)limΔt→0d0Δt+d1+d2Δt+d3Δt2+…=limΔt→0m+b1cΔt+b2kΔt2εut+2ΔtΔt=mlimΔt→0εut+2ΔtΔt=0.Equation (15) is satisfied only when(16)d0=d1=0.Combining (16) and (13) yields(17)1-b3-b6=0,2-b3=0,b1-b4-b7=0.Similarly, (8) has second-order accuracy only when d0=d1=d2=0. This means that(18)1-b3-b6=0,2-b3=0,b1-b4-b7=0,2-0.5b3-b10-b9-b11=0,2b1-b4-b9-b10-b11=0,-b9+b2-b8-b11-b5-b10=0.Solving (18) yields(19)b1=1-b7,b2=1+b8+b5,b3=2,b4=1-2b7,b6=-1,b11=1-b10-b9.This shows that USTEDM is second order accurate.Similarly, (8) has third-order accuracy only when d0=d1=d2=0. This means that(20)b1=0.5,b2=-b8,b3=2,b4=0,b5=-2b8-1,b6=-1,b7=0.5,b9=b11,b10=-2b11+1,b11=b11.Equation (8) has fourth-order accuracy when d0=d1=d2=d3=d4=0, which can be satisfied when(21)b1=0.5,b2=-b8,b3=2,b4=0,b5=-2b8-1,b6=-1,b7=0.5,b9=b11,b10=-2b11+1,-124b3-2b9-12b10+23=0,43b1-16b4-2b9-12b10=0,2b2-12b5-2b9-12b10=0.Equations (21) are contradictory equations, which shows that two-step displacement method expressed by (8) can have third-order accuracy at most, and (20) are adopted in the following derivation.


### 3.2. Stability Analysis

According to Lax theorem, (8) is convergent if and only if it is stable. The stability of (8) can be analyzed through studying a single degree-of-freedom (SDOF) system. For a SDOF system, (8) can be rewritten as(22)ut+2Δtut+Δt=b3m+b4cΔt+b5kΔt2m+b1cΔt+b2kΔt2b6m+b7cΔt+b8kΔt2m+b1cΔt+b2kΔt210ut+Δtut+b9ft+2ΔtΔt2+b10ft+ΔtΔt2+b11ft+ΔtΔt2m+b1cΔt+b2kΔt20.Amplification matrixA in (22) determining (8)’s stability is(23)A=A11A12A21A22=b3m+b4cΔt+b5kΔt2m+b1cΔt+b2kΔt2b6m+b7cΔt+b8kΔt2m+b1cΔt+b2kΔt210.Substitutingc=2ξωm,k=mω2,ωΔt=Ω(ξ-damp ratio,ω-natural frequency) into (23) yields(24)A=A11A12A21A22=b3+2b4ξΩ+b5Ω21+2b1ξΩ+b2Ω2b6+2b7ξΩ+b8Ω21+2b1ξΩ+b2Ω210.The characteristic equation of matrixA can be written as(25)det⁡A-λI=λ2-2A1λ+A2=0, where I is the identity matrix, λ denotes eigenvalue of A, and(26)A1=12tr⁡A=A11+A222=12b3+2b4ξΩ+b5Ω21+2b1ξΩ+b2Ω2,A2=A=-b6+2b7ξΩ+b8Ω21+2b1ξΩ+b2Ω2.The stability of (8) depends on spectral radius ρ(A) of matrix A and (8) is stable only when(27)ρA⩽1.Becauseρ(A) is determined by A1,A2, the stability domain of (8) can be expressed by functions A1 and A2. In similar way, the study of stability made by Hilber and Hughes [20] is performed. Firstly, the boundary line of stability domain of (8) where ρ(A)=1 is to be derived. ρ(A)=1 means that(28)λmax=1,where λmax is the maximum amplitude of λ. Complex value λ whose amplitude is 1 can be expressed by λ=eiφ, where i=-1 and φ is argument of λ. Substituting λ=eiφ into (25) yields(29)ei2φ-2A1eiφ+A2=0.The identities includingeiφ=cos⁡φ+isin⁡φ,cos⁡2φ=2cos2⁡φ-1,sin⁡2φ=2sin⁡φcos⁡φ are used and (29) can be rewritten as(30)2cos⁡φcos⁡φ-A1+A2-1+i2sin⁡φcos⁡φ-A1=0.Equation (30) is satisfied only when(31)2cos⁡φcos⁡φ-A1+A2-1=0,2sin⁡φcos⁡φ-A1=0.The boundary lines of stability domain of (8) can be determined by (31) when φ∈[0,2π]. It can be derived that when 0<φ<π,π<φ<2π, (31) are satisfied when(32)A1=cos⁡φ,A2=1.Whenφ=0,φ=2π, (31) are satisfied when(33)1-2A1+A2=0.Whenφ=π, (31) are satisfied when(34)1+2A1+A2=0.The boundary of stability domain of (8) is then made up with three lines described by (32), (33), and (34), respectively, in A1-A2 plane as shown in Figure 1. Because ρ(A) is continuous function of A1A2 and when A1=A2=0,ρ(A)=0<1, the stability domain of (8) is the inner part of the triangle as shown in Figure 1 and can be expressed by(35)-1≤A2≤1,-121+A2≤A1≤121+A2.Figure 1
Stability domain of (8).


### 3.3. Final Form of Two-Step Displacement Method with Third-Order Accuracy

Combining (20) and (8) yields(36)m+0.5cΔt-b8kΔt2ut+2Δt=2m+-2b8-1kΔt2ut+Δt+-1m+0.5cΔt+b8kΔt2ut+b11ft+2ΔtΔt2+-2b11+1ft+ΔtΔt2+b11ftΔt2.Combining (20), (26), and (35) yields, after simplification, the stability domain of (36):(37)1-b8Ω2⩾0,1+ξΩ-b8Ω2>0,4+-4b8-1Ω2⩾0.According to (36) and (37), the two-step displacement methods satisfying different requirements can be obtained. The explicitness of method can be obtained only when b11=0 and the unconditional stability can be obtained when (37) are always satisfied when Ω∈[0,∞). It can be found that (37) are always satisfied only when(38)b8⩽-0.25.Substitutingb8=-0.25,b11=0 into (36) yields(39)m+0.5cΔt+0.25kΔt2ut+2Δt=2m-0.5kΔt2ut+Δt+-m+0.5cΔt-0.25kΔt2ut+ft+ΔtΔt2.Equation (39) shows the final form of USTEDM, which is an unconditionally stable two-step explicit integration algorithm. USTEDM is third order accurate. The analysis shows that USTEDM has very similar numerical attributes as Newmark’s average acceleration method. Figure 2 shows that the spectral radius of amplification matrix of USTEDM is 1, which is the same as that of Newmark’s average acceleration method.Figure 2
Comparison of spectral radius.Figure3 shows that the numerical damping of USTEDM is zero, which is the same as that of Newmark’s average acceleration method. Figure 4 shows the relative period error of USTEDM is minimal compared with other methods including Houbolt Method, Wilson-θ Method (θ=1.4), and Generalized α Method (αm = 0.3, αf = 0.35). Though with similar numerical attributes as Newmark’s average acceleration method, USTEDM is explicit and has great computation advantage over the implicit Newmark’s average acceleration method.Figure 3
Comparison of numerical damping ratio.Figure 4
Comparison of relative period error.The applicability of the USTEDM for nonlinear problem is also investigated. Compared with other unconditional stable explicit algorithms [17–19], whose main calculation amount includes at least solving two-reverse matrix of nondiagonal matrix for displacement and velocity, respectively, USTEDM shows higher calculation efficiency since its main calculation amount includes solving only one-reverse matrix of nondiagonal matrix for displacement. In addition, USTEDM shows minimal memory requirements, since its calculation only involves displacement when only displacement-related nonlinearity occurs. This shows that USTEDM can be viewed as a complementary or improved form of conditional stable explicit CDM.


### 3.4. Implementing USTEDM for Nonlinear Problem

The procedure to use USTEDM is similar to that of CDM and is outlined as follows.(A)
Initial input and calculation:
(1)
Input the initial stiffness matrixK0, mass matrix M, and damping matrix C.(2)
Input initial displacement vectoru0 and initial velocity vector v0, solving initial acceleration vector a0:  a0=M-1(f0-K0u0-Cv0).(3)
Select time step durationΔt.(4)
Calculate(40)u-Δt=u0-v0Δt+12a0Δt2.(B)
For each time step (from time pointt to time point t+Δt):
(1)
based on knownut, ut-Δt, Kt, and outer force vector ft, calculate the “effective” load vector f-t+Δt;(41)f-t+Δt=2M-0.5KtΔt2ut+-M+0.5CΔt-0.25KtΔt2ut-Δt+ftΔt2(2)
form “effective” mass matrix:M^=M+0.5CΔt+0.25KtΔt2;(3)
calculateut+Δt=M^-1f^t+Δt;(4)
if necessary, calculate velocity vector and acceleration vector at time pointt+Δt:(42)vt+Δt=3ut+Δt-4ut+ut-Δt2Δt,at+Δt=M-1ft+Δt-Kt+Δtut+Δt-Cvt+Δt.


### 3.5. Numerical Example

In this example [7, 17], the forced vibration response of an 200-degree-of-freedom spring-mass system as shown in Figure 5 is examined using the proposed numerical method. The structural properties of this system are assumed to be mi=100 kg and ki=107[1-(ui-ui-1)2] N/m in which i=1,2,…,200. This system is subjected to a ground acceleration of 10sin⁡(πt). The lowest and highest natural frequency is 2.62 rad/s and 632.4 rad/s, respectively [17]. According to the highest natural frequency, a time step duration of Δt=0.001 s is obtained from the Newmark’s average acceleration method and is selected to be the exact solution. Displacement responses versus time of the largest degree of freedom are plotted in Figure 6. Analysis is accomplished by considering time step duration to be equal to 0.02 s. The displacement errors are presented in Figure 7. It can be seen that here USTEDM provides the most accurate displacement solutions, compared with Houbolt Method and commonly used Newmark’s average acceleration method. More importantly, without needs of iteration calculation, the unconditionally stable explicit USTEDM has much less time consumption of calculation, compared with those unconditionally stable implicit methods, which are listed in Table 1.Table 1
Total time consumption (sec).
USTEDM Houbolt Average acceleration 0.4825 353.0519 292.1621Figure 5
200-DOF spring-mass system.Figure 6
Displacement responses of 200th DOF system.Figure 7
Comparison of displacement errors of 200-DOF system.


## 3.1. Accuracy Analysis

The accuracy of USTEDM can be examined using convergence test. Equation (9) shows the Taylor series expansion of the displacement, velocity, and acceleration of the structure at time step of t+2Δt, t+Δt, respectively.(9)ut+Δt=ut+vtdt+12atdt2+16a˙tdt3+⋯,vt+Δt=vt+atdt+12a˙tdt2+16a¨tdt3+⋯,at+Δt=at+a˙tdt+12a¨tdt2+16a⃛tdt3+⋯,ut+2Δt=ut+2vtdt+2atdt2+43a˙tdt3+⋯,vt+2Δt=vt+2atdt+2a˙tdt2+43a¨tdt3+⋯,at+2Δt=at+2a˙tdt+2a¨tdt2+43a⃛tdt3+⋯.Similarly, the dynamic equilibriums att+2Δt,t+Δt,t are expressed in(10)ft+2Δt=mat+2Δt+cvt+2Δt+kut+2Δt,ft+Δt=mat+Δt+cvt+Δt+kut+Δt,ft=mat+cvt+kut.Substituting (9) into (10) yields(11)ft+2Δt=mat+2a˙tdt+⋯+cut+2vtdt+⋯+kut+2vtdt+⋯,ft+Δt=mat+a˙tdt+⋯+cvt+atdt+⋯+kut+vtdt+⋯,ft=mat+cvt+kut.Substituting (9)–(11) into (8) yields(12)m+b1cΔt+b2kΔt2εut+2Δt=m+b1cΔt+b2kΔt2ut+2Δt-b3m+b4cΔt+b5kΔt2ut+Δt-b6m+b7cΔt+b8kΔt2ut-b9ft+2ΔtΔt2-b10ft+ΔtΔt2-b11ftΔt2=d0+d1Δt+d2Δt2+d3Δt3+⋯,where(13)d0=1-b3-b6utm,d1=2-b3vtm+b1-b4-b7utc,d2=2-0.5b3-b9-b10-b11atm+2b1-b4-b9-b10-b11vtc+b2-b5-b8-b9-b10-b11utk,d3=43-16b3-2b9-b10a˙tm+2b1-0.5b4-2b9-b10atc+2b2-b5-2b9-b10vtk,d4=23-124b3-2b9-12b10a¨tm+43b1-16b4-2b9-12b10a˙tc+2b2-12b5-2b9-12b10atk,d5=415-1120b3-43b9-16b10a⃛tm+23b1-124b4-43b9-16b10a¨tc+43b2-16b5-43b9-16b10a˙tk⋯.USTEDM having first-order accuracy means(14)limΔt→0εut+2ΔtΔt=0.Combining (12) and (14) leads to(15)limΔt→0d0Δt+d1+d2Δt+d3Δt2+…=limΔt→0m+b1cΔt+b2kΔt2εut+2ΔtΔt=mlimΔt→0εut+2ΔtΔt=0.Equation (15) is satisfied only when(16)d0=d1=0.Combining (16) and (13) yields(17)1-b3-b6=0,2-b3=0,b1-b4-b7=0.Similarly, (8) has second-order accuracy only when d0=d1=d2=0. This means that(18)1-b3-b6=0,2-b3=0,b1-b4-b7=0,2-0.5b3-b10-b9-b11=0,2b1-b4-b9-b10-b11=0,-b9+b2-b8-b11-b5-b10=0.Solving (18) yields(19)b1=1-b7,b2=1+b8+b5,b3=2,b4=1-2b7,b6=-1,b11=1-b10-b9.This shows that USTEDM is second order accurate.Similarly, (8) has third-order accuracy only when d0=d1=d2=0. This means that(20)b1=0.5,b2=-b8,b3=2,b4=0,b5=-2b8-1,b6=-1,b7=0.5,b9=b11,b10=-2b11+1,b11=b11.Equation (8) has fourth-order accuracy when d0=d1=d2=d3=d4=0, which can be satisfied when(21)b1=0.5,b2=-b8,b3=2,b4=0,b5=-2b8-1,b6=-1,b7=0.5,b9=b11,b10=-2b11+1,-124b3-2b9-12b10+23=0,43b1-16b4-2b9-12b10=0,2b2-12b5-2b9-12b10=0.Equations (21) are contradictory equations, which shows that two-step displacement method expressed by (8) can have third-order accuracy at most, and (20) are adopted in the following derivation.


## 3.2. Stability Analysis

According to Lax theorem, (8) is convergent if and only if it is stable. The stability of (8) can be analyzed through studying a single degree-of-freedom (SDOF) system. For a SDOF system, (8) can be rewritten as(22)ut+2Δtut+Δt=b3m+b4cΔt+b5kΔt2m+b1cΔt+b2kΔt2b6m+b7cΔt+b8kΔt2m+b1cΔt+b2kΔt210ut+Δtut+b9ft+2ΔtΔt2+b10ft+ΔtΔt2+b11ft+ΔtΔt2m+b1cΔt+b2kΔt20.Amplification matrixA in (22) determining (8)’s stability is(23)A=A11A12A21A22=b3m+b4cΔt+b5kΔt2m+b1cΔt+b2kΔt2b6m+b7cΔt+b8kΔt2m+b1cΔt+b2kΔt210.Substitutingc=2ξωm,k=mω2,ωΔt=Ω(ξ-damp ratio,ω-natural frequency) into (23) yields(24)A=A11A12A21A22=b3+2b4ξΩ+b5Ω21+2b1ξΩ+b2Ω2b6+2b7ξΩ+b8Ω21+2b1ξΩ+b2Ω210.The characteristic equation of matrixA can be written as(25)det⁡A-λI=λ2-2A1λ+A2=0, where I is the identity matrix, λ denotes eigenvalue of A, and(26)A1=12tr⁡A=A11+A222=12b3+2b4ξΩ+b5Ω21+2b1ξΩ+b2Ω2,A2=A=-b6+2b7ξΩ+b8Ω21+2b1ξΩ+b2Ω2.The stability of (8) depends on spectral radius ρ(A) of matrix A and (8) is stable only when(27)ρA⩽1.Becauseρ(A) is determined by A1,A2, the stability domain of (8) can be expressed by functions A1 and A2. In similar way, the study of stability made by Hilber and Hughes [20] is performed. Firstly, the boundary line of stability domain of (8) where ρ(A)=1 is to be derived. ρ(A)=1 means that(28)λmax=1,where λmax is the maximum amplitude of λ. Complex value λ whose amplitude is 1 can be expressed by λ=eiφ, where i=-1 and φ is argument of λ. Substituting λ=eiφ into (25) yields(29)ei2φ-2A1eiφ+A2=0.The identities includingeiφ=cos⁡φ+isin⁡φ,cos⁡2φ=2cos2⁡φ-1,sin⁡2φ=2sin⁡φcos⁡φ are used and (29) can be rewritten as(30)2cos⁡φcos⁡φ-A1+A2-1+i2sin⁡φcos⁡φ-A1=0.Equation (30) is satisfied only when(31)2cos⁡φcos⁡φ-A1+A2-1=0,2sin⁡φcos⁡φ-A1=0.The boundary lines of stability domain of (8) can be determined by (31) when φ∈[0,2π]. It can be derived that when 0<φ<π,π<φ<2π, (31) are satisfied when(32)A1=cos⁡φ,A2=1.Whenφ=0,φ=2π, (31) are satisfied when(33)1-2A1+A2=0.Whenφ=π, (31) are satisfied when(34)1+2A1+A2=0.The boundary of stability domain of (8) is then made up with three lines described by (32), (33), and (34), respectively, in A1-A2 plane as shown in Figure 1. Because ρ(A) is continuous function of A1A2 and when A1=A2=0,ρ(A)=0<1, the stability domain of (8) is the inner part of the triangle as shown in Figure 1 and can be expressed by(35)-1≤A2≤1,-121+A2≤A1≤121+A2.Figure 1
Stability domain of (8).


## 3.3. Final Form of Two-Step Displacement Method with Third-Order Accuracy

Combining (20) and (8) yields(36)m+0.5cΔt-b8kΔt2ut+2Δt=2m+-2b8-1kΔt2ut+Δt+-1m+0.5cΔt+b8kΔt2ut+b11ft+2ΔtΔt2+-2b11+1ft+ΔtΔt2+b11ftΔt2.Combining (20), (26), and (35) yields, after simplification, the stability domain of (36):(37)1-b8Ω2⩾0,1+ξΩ-b8Ω2>0,4+-4b8-1Ω2⩾0.According to (36) and (37), the two-step displacement methods satisfying different requirements can be obtained. The explicitness of method can be obtained only when b11=0 and the unconditional stability can be obtained when (37) are always satisfied when Ω∈[0,∞). It can be found that (37) are always satisfied only when(38)b8⩽-0.25.Substitutingb8=-0.25,b11=0 into (36) yields(39)m+0.5cΔt+0.25kΔt2ut+2Δt=2m-0.5kΔt2ut+Δt+-m+0.5cΔt-0.25kΔt2ut+ft+ΔtΔt2.Equation (39) shows the final form of USTEDM, which is an unconditionally stable two-step explicit integration algorithm. USTEDM is third order accurate. The analysis shows that USTEDM has very similar numerical attributes as Newmark’s average acceleration method. Figure 2 shows that the spectral radius of amplification matrix of USTEDM is 1, which is the same as that of Newmark’s average acceleration method.Figure 2
Comparison of spectral radius.Figure3 shows that the numerical damping of USTEDM is zero, which is the same as that of Newmark’s average acceleration method. Figure 4 shows the relative period error of USTEDM is minimal compared with other methods including Houbolt Method, Wilson-θ Method (θ=1.4), and Generalized α Method (αm = 0.3, αf = 0.35). Though with similar numerical attributes as Newmark’s average acceleration method, USTEDM is explicit and has great computation advantage over the implicit Newmark’s average acceleration method.Figure 3
Comparison of numerical damping ratio.Figure 4
Comparison of relative period error.The applicability of the USTEDM for nonlinear problem is also investigated. Compared with other unconditional stable explicit algorithms [17–19], whose main calculation amount includes at least solving two-reverse matrix of nondiagonal matrix for displacement and velocity, respectively, USTEDM shows higher calculation efficiency since its main calculation amount includes solving only one-reverse matrix of nondiagonal matrix for displacement. In addition, USTEDM shows minimal memory requirements, since its calculation only involves displacement when only displacement-related nonlinearity occurs. This shows that USTEDM can be viewed as a complementary or improved form of conditional stable explicit CDM.


## 3.4. Implementing USTEDM for Nonlinear Problem

The procedure to use USTEDM is similar to that of CDM and is outlined as follows.(A)
Initial input and calculation:
(1)
Input the initial stiffness matrixK0, mass matrix M, and damping matrix C.(2)
Input initial displacement vectoru0 and initial velocity vector v0, solving initial acceleration vector a0:  a0=M-1(f0-K0u0-Cv0).(3)
Select time step durationΔt.(4)
Calculate(40)u-Δt=u0-v0Δt+12a0Δt2.(B)
For each time step (from time pointt to time point t+Δt):
(1)
based on knownut, ut-Δt, Kt, and outer force vector ft, calculate the “effective” load vector f-t+Δt;(41)f-t+Δt=2M-0.5KtΔt2ut+-M+0.5CΔt-0.25KtΔt2ut-Δt+ftΔt2(2)
form “effective” mass matrix:M^=M+0.5CΔt+0.25KtΔt2;(3)
calculateut+Δt=M^-1f^t+Δt;(4)
if necessary, calculate velocity vector and acceleration vector at time pointt+Δt:(42)vt+Δt=3ut+Δt-4ut+ut-Δt2Δt,at+Δt=M-1ft+Δt-Kt+Δtut+Δt-Cvt+Δt.


## 3.5. Numerical Example

In this example [7, 17], the forced vibration response of an 200-degree-of-freedom spring-mass system as shown in Figure 5 is examined using the proposed numerical method. The structural properties of this system are assumed to be mi=100 kg and ki=107[1-(ui-ui-1)2] N/m in which i=1,2,…,200. This system is subjected to a ground acceleration of 10sin⁡(πt). The lowest and highest natural frequency is 2.62 rad/s and 632.4 rad/s, respectively [17]. According to the highest natural frequency, a time step duration of Δt=0.001 s is obtained from the Newmark’s average acceleration method and is selected to be the exact solution. Displacement responses versus time of the largest degree of freedom are plotted in Figure 6. Analysis is accomplished by considering time step duration to be equal to 0.02 s. The displacement errors are presented in Figure 7. It can be seen that here USTEDM provides the most accurate displacement solutions, compared with Houbolt Method and commonly used Newmark’s average acceleration method. More importantly, without needs of iteration calculation, the unconditionally stable explicit USTEDM has much less time consumption of calculation, compared with those unconditionally stable implicit methods, which are listed in Table 1.Table 1
Total time consumption (sec).
USTEDM Houbolt Average acceleration 0.4825 353.0519 292.1621Figure 5
200-DOF spring-mass system.Figure 6
Displacement responses of 200th DOF system.Figure 7
Comparison of displacement errors of 200-DOF system.


## 4. Conclusions

A novel two-step unconditionally stable explicit displacement method, named Unconditional Stable Two-Step Explicit Displacement Method (USTEDM), is proposed in this paper. USTEDM has third-order accuracy and is unconditionally stable. The conclusions mainly include the following:(1)
The novel way of constructing direct time integration methods from the start point of dimensional analysis is applicable, and it is believed that it will contribute to improvement of numerical calculation in structural dynamics.(2)
With no prerequisite difference assumption of velocity and acceleration, the derivation of USTEDM is completely different from any documented one. With unconditional stability, third-order accuracy, higher calculation efficiency, and lower storage requirement, USTEDM is promising in solving vibrations of structures, especially those with displacement-dependent nonlinearity.The zero numerical damping of USTEDM is its shortcoming, and the introduction of numerical damping will be studied.



---

*Source: 2831206-2016-12-08.xml*