# Two-Step Unconditionally Stable Noniterative Dissipative Displacement Method for Analysis of Nonlinear Structural Dynamics Problems

Authors: Changqing Li; Baoyi Sheng; Zhipeng Lai; Lizhong Jiang; Ping Xiang
Journal: Shock and Vibration
 (2021)

---

## Abstract

When solving structural dynamic problems, the displacement algorithm needs only calculating and storing structure’s displacements in the main calculation process, which makes the displacement algorithm have advantages over multivariable algorithms in calculation efficiency and storage requirements. By using a novel approach based on dimensional analysis firstly given by the first author, a one-parameter family of two-step unconditionally stable noniterative displacement algorithms, referred to as the CQ-2x method, is developed. Compared with other unconditionally stable noniterative multivariable algorithms such as the representative KR-α method, the proposed method has advantages in several aspects. The CQ-2x method is unconditionally stable regardless of stiffness hardening or stiffness weakening, while the KR-α method is only conditionally stable in case of stiffness hardening. The CQ-2x method needs only one solver within one time step, while the KR-α method needs two solvers within one time step, which makes the CQ-2x method show higher efficiency. Numerical examples are presented to demonstrate the potential of the proposed method.


---

## Body

## 1. Introduction

Direct integration methods are commonly used to solve structural dynamics equations. According to whether iteration is needed to solve nonlinear dynamic problems, these methods can be divided into two categories: implicit and noniterative. If the solution for the next time step depends on the responses in future time steps, then the algorithm is classified as an implicit algorithm [1, 2]. The well-known Houbolt [3], Newmark [4], Wislon [5], HHT-α [6], WBZ α [7], generalized-α [8], and Bath [9–11] algorithms are examples of implicit algorithms. When the implicit methods are used to solve nonlinear structural dynamics equations, an iterative scheme (for example, the Newton–Raphson algorithm) is required in each time step which is often complicated and time consuming.By contrast, equilibrium in noniterative algorithms is satisfied in the previous time step rather than the current step; consequently, such algorithms do not require the formation of the tangent stiffness matrix or any iterative calculations and are computationally more efficient than implicit algorithms in the nonlinear case when using same time step size.The central difference method (CDM) [1, 2] is the most common and effective noniterative method since it not only needs no nonlinear iteration but also needs no solver for nonlinear problems. However, the CDM is only conditionally stable, which means extraordinary small time step size must be used in order to make the calculation stable. The conditionally stable CDM is preferred when the time step size required for accuracy is of the same order as the stability limit of the algorithm, e.g., in wave propagation problems [2, 12]. Nevertheless, unconditionally stable noniterative algorithms are more useful for solving nonlinear structural dynamics problems.Chang’s algorithms [13–16], the CR [17] method, and the KR-α [18] method are examples of unconditionally stable noniterative algorithms. The first unconditionally stable noniterative structure-dependent integration method was developed by Chang [13]. An assessment of such structure-dependent unconditionally stable noniterative and semiexplicit (SE) algorithms for application to dynamic problems has been presented by Kolay and Ricles [19]. These authors found that an SE algorithm may generate large damping forces when a realistic time step size is used for the nonlinear transient analysis of realistic structural systems subjected to seismic excitations, whereas the KR-α method [18] may produce high-frequency overshoot. Furthermore, unconditionally stable noniterative model-dependent algorithms may become conditionally stable when severe stiffness hardening occurs, and most of them require the solution of linear equations whose coefficient matrix is nondiagonal twice within the one time step, which lowers their computational efficiency. When the number of DOFs of the structure is large, the majority of the computation time of an algorithm is consumed in solving linear equations whose coefficient matrices are nondiagonal. From the perspective of computational efficiency, it is preferable for an algorithm to solve a system of linear equations only once per time step.In conjunction with the total energy framework and a generalized time-weighted residual approach, Har and Tamma [20, 21] have provided good guidance for developing new time integration algorithms and procedures for evaluating new and existing algorithms. However, the unconditionally stable noniterative structure-dependent integration methods can hardly be evaluated by those guidances. The algorithms including the structure-dependent algorithm may be easier to be designed by using a novel approach developed by the first author of the current study based on dimensional analysis. Most of the commonly used algorithms for structural dynamics analysis consist of one motion equation and two difference equations of displacement and velocity. By contrast, the novel dimensional analysis approach presented here relies on only a simple objective function or functions (the case of multiple functions will not be discussed in this paper), and the assumptions concerning displacement, velocity, and acceleration are implicitly included in their Taylor series expansions and the dynamic equilibrium equation.When solving structural dynamic problems, according to the types of physical quantities that must be involved in the main calculation steps, the algorithm can be divided into the single variable method and multivariable method. As a single variable method, the displacement method needs only calculating and storing displacements in the calculation process, which makes the displacement method have advantages over other multivariable algorithms in calculation efficiency and storage requirements. So the novel approach was applied firstly to design the displacement method. By using this novel approach, a two-step unconditionally stable noniterative displacement algorithm (NDA) called the CQ-2x method (a two-step algorithm named after the first author) has been proposed [22]. The CQ-2x method combines unconditional stability with noniterative formulations of both displacement and velocity, the use of one solver per time step, and no overshoot in either displacement or velocity. However, CQ-2x does not offer controllable numerical dissipation. To introduce controllable numerical damping, while simultaneously preserving the previously mentioned desirable numerical characteristics of the CQ-2x algorithm, a three-step displacement algorithm called the CQ-3 method [23] with controllable numerical dissipation has also been developed. However, the need for three steps is a disadvantage that may lower the efficiency of the algorithm. Thus, the development of a family of two-step unconditionally stable NDAs with controllable numerical damping is the main task of this paper.The article begins with an objective function for a two-step NDA. Then, the preliminary form of the two-step NDA is obtained through dimensional analysis. Since there are 9 coefficients to be determined in the preliminary form of the NDA, equations for determining those coefficients from the algorithm’s consistency and stability are provided. Since the algorithm’s consistency and stability requirements cannot be used to determine all coefficients, the remaining free coefficientsb5, b7, and b8 are then determined from other characteristics of the algorithm, such as the relative period error and numerical damping ratio. A self-starting procedure is presented, and the overshoots in displacement and velocity are analyzed; subsequently, three numerical examples are provided to demonstrate the accuracy, stability, and efficiency of the proposed method.


## 2. Preliminary Form of NDA

For a single-DOF system, the linear structural dynamics equation can be written as(1)mat+cvt+kut=ft.In the equation above,m,c,andk represent the mass, damping, and stiffness, respectively; a, v, and u represent the acceleration, velocity, and displacement of the mass, respectively; f is the external force; and t is the time.When no nonlinear damping is considered, the nonlinear structural dynamics equation can be written as follows [18]:(2)mat+cvt+finu,t=fextt,where fextt is the external force.The nonlinearity of a structure may be complex. In some cases, the structural nonlinearity can be expressed in terms of nonlinear damping and nonlinear stiffness. In these cases, the nonlinear structural dynamics equation can be rewritten as(3)mat+ctvt+ktut=ft,where the nonlinear damping ct and the nonlinear stiffness kt are given by(4)kt=kut,vt,t,ct=cut,vt,t.In order to obtain the numerical solution at the discrete time point, according to the definition of the “noniterative” displacement algorithm given above, the objective of a two-step NDA is to findut+2Δt when ut+Δt, ut, ft+Δt, m, ct+Δt, kt+Δt, and Δt are given; this can be expressed by the objective function given in equation (5) where ut+2Δt is regarded as a dependent variable and ut+Δt, ut, ft+Δt, m, ct+Δt, kt+Δt, and Δt are regarded as independent variables:(5)ut+2Δt=ut+2Δtm,ct+Δt,kt+Δt,ft+Δt,ut+Δt,ut,Δt,where Δt is the time step size and the subscripts denote discrete time points. Because the displacement algorithm is to be designed, there are no velocity and acceleration occurring in equation (5). A displacement algorithm designed based on dimensional analysis to solve the objective function given in equation (5) can be written as follows:(6)m+b1ct+ΔtΔt+b2kt+ΔtΔt2ut+2Δt=b3m+b4ct+ΔtΔt+b5kt+ΔtΔt2ut+Δt+b6m+b7ct+ΔtΔt+b8kt+ΔtΔt2ut+b9ft+ΔtΔt2,where bi are undetermined coefficients. The detailed derivation of equation (6) based on dimensional analysis is presented in Appendix A. bi can be determined from the NDA’s consistency and stability requirements.


## 3. Consistency Analysis to Determinebi

The consistency of an NDA of the form given in equation (6) can be analyzed by calculating the local algorithmic error. To calculate the local algorithmic error εu2+3Δt between the exact displacement ut+2Δt and the computed displacement ut+2Δt, it is assumed that, in equation (6),(7a)ut=ut,(7b)ut+Δt=ut+Δt,(7c)ut+2Δt=ut+2Δt−εut+2Δt,(7d)ct+Δt=ct+Δt,(7e)kt+Δt=kt+Δt,(7f)ft+Δt=ft+Δt,where a quantity with brackets “()” denotes an exact quantity and a quantity without brackets “()” but with a subscript denotes a computed quantity.Substituting equations (7a)–(7f) into equation (6) yields the following:(8)m+b1ct+ΔtΔt+b2kt+ΔtΔt2εut+2Δt=m+b1ct+ΔtΔt+b2kt+ΔtΔt2ut+2Δt−b3m+b4ct+ΔtΔt+b5kt+ΔtΔt2ut+Δt−b6m+b7ct+ΔtΔt+b8kt+ΔtΔt2ut−b9ft+ΔtΔt2.Taylor series expansions ofut+Δt, ut+2Δt, vt+Δt, and at+Δt can be written as follows:(9a)ut+Δt=ut+vtΔt+12atΔt2+16a˙tΔt3+⋯,(9b)ut+2Δt=ut+2vtΔt+2atΔt2+43a˙tΔt3+⋯,(9c)vt+Δt=vt+atΔt+12a˙tΔt2+16a¨tΔt3+⋯,(9d)at+Δt=at+a˙tΔt+12a¨tΔt2+16a⃛tΔt3+⋯,where a superscripted dot denotes differentiation with respect to time. The nonlinear dynamic equilibrium of equation (3) in step t+Δt is written as(10)ft+Δt=mat+Δt+ct+Δtvt+Δt+kt+Δtut+Δt.Substituting equations (9a), (9c), and (9d) into equation (10) yields(11)ft+Δt=mat+a˙tΔt+⋯+ct+Δtvt+atΔt+⋯+kt+Δtut+vtΔt+⋯.Substituting equations (9a), (9b), and (11) into equation (8) yields(12)m+b1ct+ΔtΔt+b2kt+ΔtΔt2εut+2Δt=e0+e1Δt+e2Δt2+e3Δt3+e4Δt4+e5Δt5⋯,where the expressions for e0, e1, e2, e3, e4, and e5 are given in Appendix B.Rewriting equation (12) yields(13)limΔt⟶0εut+2Δt=limΔt⟶0e0+e1Δt+e2Δt2+e3Δt3+e4Δt4+e5Δt5⋯m+b1ct+ΔtΔt+b2kt+ΔtΔt2,=limΔt⟶0e0+e1Δt+e2Δt2+e3Δt3+e4Δt4+e5Δt5⋯m.The statement that the NDA has a nominalnth order time accuracy of displacement means that the following holds:(14)limΔt⟶0εut+2ΔtΔtn=0.The combination of equations (13) and (14) leads to(15)limΔt⟶0εut+2ΔtΔtn=1me0Δtn+e1Δtn−1+e2Δtn−2+e3Δtn−3+⋯+=0.Substitutingn=2 into equation (15) yields the condition that must be met for the NDA to achieve a nominal second-order time accuracy of displacement:(16)limΔt⟶0εut+2ΔtΔt2=1me0Δt2+e1Δt+e2+e3Δt+⋯+=0.Equation (16) is satisfied only when e0=e1=e2=0, and the conditions for which can be obtained by simultaneously solving equations (B.2a), (B.2b), and (B.2c) (the detailed derivation can be seen in Appendix B):(17)b1=1−b7,b2=1+b8+b5,b3=2,b4=1−2b7,b6=−1,b9=1.Equation (17) ensures that the NDA will achieve a nominal second-order time accuracy of displacement.To confirm the nominal second-order time accuracy of displacement achieved with equation (17), we can substitute equation (17) into equation (12) and obtain the following equation:(18)εut+2Δt=−Δt31201+b5+2b8ωt+Δt2vt+1201−2b7atξt+Δtωt+Δt+ϑΔt4120b5+b8+1Ωt+Δt2+2401−b7ξt+ΔtΩt+Δt+120,where(19)ωt+Δt=kt+Δtm,ξt+Δt=ct+Δt2mωt+Δt,Ωt+Δt=ωt+ΔtΔt.The substitution of equation (17) into equation (6) yields an NDA with a nominal second-order time accuracy of displacement:(20)m+1−b7ct+ΔtΔt+1+b8+b5kt+ΔtΔt2ut+2Δt=2m+1−2b7ct+ΔtΔt+b5kt+ΔtΔt2ut+Δt+−m+b7ct+ΔtΔt+b8kt+ΔtΔt2ut+ft+ΔtΔt2.It can be seen from the derivation of equation (20) that the NDA is derived from only one objective function, equation (5). The possible assumed relationships among the displacement, velocity, and acceleration are all implied in the Taylor series expansions (as seen in equations (9a)–(9d)) and the nonlinear dynamic equilibrium equation (as seen in equation (10)). Newmark’s assumptions (with either constant coefficients or structure-dependent coefficients) concerning displacement and velocity are not needed for the NDA family derived in this paper although they are needed for most other algorithms; thus, this derivation is more comprehensive and allows no promising algorithms to be missed when the objective function is determined. It should be noted here that the NDA accuracy analysis applies to the case of nonlinear structural dynamics since nonlinear equation (10) is used in the derivation instead of linear equation (1), which is also a point of difference with respect to other algorithms.When calculation of the velocity is necessary, it can be achieved with acceptable accuracy by means of any applicable displacement difference approximation. For example, one can use, but is not limited to, the displacement difference approximation used in the central-eccentric difference method [24]:(21)vt+2Δt=3ut+2Δt−4ut+Δt+ut2Δt.Equation (21) has a nominal time accuracy of velocity of at least first order, and the local algorithmic velocity error εvt+2Δt=vt+2Δt−vt+2Δt when this approximation is applied in the proposed CQ-2x method is(22)εvt+2Δt=−Δt280a˙t+3601+b5+2b8ωt+Δt2vt+3601−2b7ξt+Δtωt+Δtat+φΔt3240b5+b8+Ω2+4801−b7ξΩ+240.After the noniterative calculation ofut+2Δt and vt+2Δt, at+2Δt can be calculated from the nonlinear dynamic equilibrium (as given in equation (3)) at step t+2Δt:(23)at+2Δt=1mft+2Δt−ct+2Δtvt+2Δt−kt+2Δtut+2Δt.Equation (23) has a nominal time accuracy of acceleration of at least first order, and when this approximation is applied in the CQ-2x method, the local algorithmic acceleration error εat+2Δt=at+2Δt−at+2Δt can be written as(24)εat+2Δt=Δt2ξt+2Δtωt+2Δt2a˙t+91+b5+2b8ωt+Δt2vt+91−2b7ξt+Δtωt+Δtat+ϑΔt331+b5+b8Ωt+Δt2+61−b7ξt+ΔtΩt+Δt+3,where(25)ωt+2Δt=kt+2Δtm,ξt+2Δt=ct+2Δt2mωt+2Δt.The accuracy expressed above is called the nominal accuracy because the derivation of the displacement consistency assumes that the values ofu at time steps t and t+Δt are known exactly, which cannot be ensured, meaning that this assumption leads to additional errors. However, it can be concluded that an NDA described by equations (20), (21), and (23) has a global time accuracy of at least first order, as will be validated by Example 1.


## 4. Stability Analysis to Determinebi

An algorithm is convergent if and only if it is consistent and stable. The consistency of the proposed NDA family has been analyzed above. In this section, the stability of an NDA with the form given in equation (20) will be studied. It should again be noted here that the NDA consistency analysis applies to the solution of nonlinear structural dynamics problems of the form given in equation (3) because equation (3), not equation (1), is used throughout the entire derivation. For the same reason, the following NDA stability analysis also applies to the solution of nonlinear structural dynamics problems of the form given in equation (3). The NDA expressed in equation (20) can be rewritten as follows:(26)ut+2Δtut+Δt=Aut+Δtut+b9ft+ΔtΔt2m+1−b7ct+ΔtΔt+1+b8+b5kt+ΔtΔt20,where the amplification matrix A is(27)A=A11A1210,(28a)%A11=2m+1−2b7ct+ΔtΔt+b5kt+ΔtΔt2m+1−b7ct+ΔtΔt+1+b8+b5kt+ΔtΔt2,(28b)A12=−m+b7ct+ΔtΔt+b8kt+ΔtΔt2m+1−b7ct+ΔtΔt+1+b8+b5kt+ΔtΔt2.A stability study similar to the one conducted by Hilber and Hughes [25] is performed (the detailed derivation can be seen in Appendix C). The stability conditions for a first-order time-accurate NDA can be written as follows:(29a)−1≤A2≤1,(29b)−121+A2≤A1≤121+A2,where(30a)A1=12trA=A11+A222=122m+1−2b7ct+ΔtΔt+b5kt+ΔtΔt2m+1−b7ct+ΔtΔt+1+b8+b5kt+ΔtΔt2,(30b)A2=A=−−m+b7ct+ΔtΔt+b8kt+ΔtΔt2m+1−b7ct+ΔtΔt+1+b8+b5kt+ΔtΔt2.By consideringωt+Δt=kt+Δtm−1, ct+Δt=2mωt+Δtξt+Δt, and ωt+ΔtΔt=Ωt+Δt, where ωt+Δt is the natural frequency and ξt+Δt is the damping ratio; equations (30a) and (30b) at time t+Δt can be rewritten as follows:(31a)A1=12trA=A11+A222=122+21−2b7ξt+ΔtΩt+Δt+b5Ωt+Δt21+21−b7ξt+ΔtΩt+Δt+1+b8+b5Ωt+Δt2.(31b)A2=A=−−1+2b7ξt+ΔtΩt+Δt+b8Ωt+Δt21+21−b7ξt+ΔtΩt+Δt+1+b8+b5Ωt+Δt2.Substituting equations (31a) and (31b) into equations (29a) and (29b) yields the following:(32a)1+1−b72ξt+ΔtΩt+Δt+1+b8+b5Ωt+Δt2≥0,(32b)4+2−4b72ξt+ΔtΩt+Δt+1+2b5Ωt+Δt2≥0,(32c)2ξt+ΔtΩt+Δt+1+2b8+b5Ωt+Δt2≥0.Equations (32a)–(32c) give the stability conditions for the NDA that apply when solving a nonlinear structural dynamics problem of the form given in equation (3) with at least first-order time accuracy. Equation (20) shows that the CDM is only one specific kind of NDA with b5=−1, b7=1/2, and b8=0. To check the correctness of the NDA stability conditions given in equations (32a)–(32c), one can solve for the domain of stability of the CDM simply by substituting b5=−1, b7=1/2, and b8=0 into equations (32a)–(32c). Then, equations (32a)–(32c) take the following form:(33a)1+ξt+ΔtΩt+Δt≥0,(33b)4−Ωt+Δt2≥0,(33c)2ξt+ΔtΩt+Δt≥0.It can be seen that equations (33a)–(33c) yield the correct stability domain for the CDM, Ωt+Δt2≤4. Since the stiffness ratio δ (where δ<1 indicates stiffness softening and δ>1 indicates stiffness hardening) does not appear in the amplification matrix A of CQ-2x but does appear in the amplification matrices of other algorithms (for example, Chang’s algorithms [16]), which causes them to become unstable when stiffness hardening occurs, it may be concluded that when used to solve the nonlinear structural dynamics equation given in equation (3), the CQ-2x method is unconditionally stable regardless of whether stiffness hardening or softening occurs as long as equations (32a)–(32c) are always satisfied when Ω∈0,∞ and ξ∈0,∞, which will be validated by Example 2.


## 5. CQ-2x Method with Controllable Numerical Dissipation

Up to this point, three coefficients have remained undetermined:b5, b7, and b8. These coefficients can be determined from the algorithm’s characteristics of unconditional stability and acceptable numerical dispersion and dissipation.


### 5.1. Determination ofb7

The unconditional stability of an NDA of the form given in equation (20) means that when Ω∈0,+∞ and ξ∈0,+∞, equations (32a)–(32c) are always satisfied, which requires that the following must hold:(34)1−b7≥0,1+b8+b5≥0,2−4b7≥0,1+2b5≥0,1+2b8+b5≥0.An intermediate variablex is introduced:(35)x=1+2b8+b5.Consequently,b8 can be written as(36)b8=12x−b5−1.Then, simplifying equation (34) leads to(37)b7≤0.5,b8≤x,b5≥−0.5,x≥0.An analysis shows that, for a damping structure, a largerb7 leads to more satisfactory algorithmic characteristics; thus, b7 is assigned its largest possible value of 0.5 according to equation (37):(38)b7=0.5.Substituting equations (36) and (38) into equation (20) yields the preliminary algorithmic form of the CQ-2x method:(39)m+0.5ct+ΔtΔt+12x+b5+1kt+ΔtΔt2ut+2Δt=2m+b5kt+ΔtΔt2ut+Δt+−m+0.5ct+ΔtΔt+12x−b5−1kt+ΔtΔt2ut+ft+ΔtΔt2.It can be seen from the above analysis that there are only two coefficients,b5 and x, that are undetermined in the preliminary form of the CQ-2x method. These coefficients will be determined by introducing controllable numerical dissipation into CQ-2x.


### 5.2. Determination ofb5 and x

Whenξ=0, according to equation (36), equations (31a) and (31b) can be written as(40a)%A1=A11+A222=1+1/2b5Ω21+1/2b5+1/2x+1/2Ω2,(40b)A2=1+1/2+1/2b5−1/2xΩ21+1/2+1/2b5+1/2xΩ2.The effects ofb5 and x on the relative period error of the preliminary CQ-2x method have been analyzed. The results when Δt/T=0.1, b5=−0.5,−0.3,0,0.3,1,2,and4, and x∈0,4 are plotted in Figure 2. It can be seen that, as b5 and x increase, the relative period error also increases. Thus, to achieve the best accuracy, x and b5 should be assigned their lowest possible values of 0 and −0.5, respectively. The analysis shows that when x=0 and b5=−0.5, the preliminary CQ-2x method has a nominal third-order time accuracy of displacement and a global second-order time accuracy; this will be validated by Example 1. However, if numerical dissipation is needed to filter out spurious vibrations in high-vibration modes, neither x=0 nor b5=−0.5 is appropriate.Figure 1
NDA stability domain.Figure 2
Variations in the relative period error of the preliminary CQ-2x method with b5 and x when Δt/T=0.1.It can be proven that whenx=0, the numerical damping ratio is also ξ¯=0, as shown in the following equation:(41)x=0,A2≥A12,A2=1+1/2+1/2b5Ω21+1/2+1/2b5Ω2=1,ρA=A2=1,ξ¯=0.It can also be proven that whenb5=−0.5, the numerical damping ratio ξ¯ of very high vibration modes Ω⟶∞ is also zero, as shown in the following equation:(42)b5=−0.5,A2≤A12,limΩ⟶∞ρA=limΩ⟶∞A1−A12−A2=−1/4−1/2x1/2x+1/4=1,limΩ⟶∞ξ¯=0.To introduce an appropriate level of numerical dissipation,b5 and x must be greater than their lowest values of 0 and −0.5, respectively. The numerical damping ratio depends on the spectral radius of the amplification matrix. It may be preferable for the spectral radius of the amplification matrix A of the preliminary CQ-2x method to decrease with increasing Ω with no bifurcation. According to equations (D1a) and (D1b), no bifurcation occurs when A2≥A12, which requires the following:(43)A2−A12=−x24−b52−14Ω4+Ω2≥0,whenΩ∈0,+∞.Equation (43) is satisfied when(44)b5≥x22−12.Since a lower value ofb5 results in a higher accuracy, as mentioned before, b5 is assigned the lowest possible value according to equation (44):(45)b5=x22−12.Substituting equation (45) into equation (39) yields(46)m+12ct+ΔtΔt+x+124kt+ΔtΔt2ut+2Δt=2m+x2−12kt+ΔtΔt2ut+Δt+−m+12ct+ΔtΔt−x−124kt+ΔtΔt2ut+ft+ΔtΔt2.Equation (46) is the final algorithmic form of displacement calculation in the CQ-2x method.To determinex, while simultaneously confirming the unconditional stability of the CQ-2x method, consider the following two extreme cases: Ω⟶0 and Ω⟶∞. In accordance with equations (31a) and (31b), (36), (45), and (46) and considering the limits, it can be shown that equation (C1) reduces to the following:(47a)λ−12=0,whenΩ⟶0,(47b)λ−x−1x+12=0,whenΩ⟶∞.It can be calculated from equation (47b) that(48)x=1−ρ∞1+ρ∞.Figures3 and 4 show how the spectral radius of the CQ-2x method varies with ρ∞ and Ω; the findings are fully consistent with the above analysis.Figure 3
Spectral radius of CQ-2x when ξ=0.Figure 4
Spectral radius of CQ-2x when ξ=0.05.


### 5.3. The Proposed Algorithm

Now that all coefficients in the preliminary NDA form presented in equation (6) have been determined, a family of unconditionally stable noniterative algorithms, collectively referred to as the CQ-2x method with one-parameter x, can be obtained by combining equations (21), (23), (46), and (48), which are rewritten as(49a)m+12ct+ΔtΔt+x+124kt+ΔtΔt2ut+2Δt=2m+x2−12kt+ΔtΔt2ut+Δt+−m+12ct+ΔtΔt−x−124kt+ΔtΔt2ut+ft+ΔtΔt2,(49b)vt+2Δt=3ut+2Δt−4ut+Δt+ut2Δt,(49c)at+2Δt=1mft+2Δt−ct+2Δtvt+2Δt−kt+2Δtut+2Δt,where(49d)x=1−ρ∞1+ρ∞.


### 5.4. Stability for Systems with Nonlinear Stiffness

It has been demonstrated that the CQ-2x method is unconditionally stable when it is used to solve the nonlinear structural dynamics equation given in equation (3). Here, the stability of CQ-2x when it is used to solve the nonlinear structural dynamics equation given in equation (2) will be analyzed. When a general nonlinear structural dynamics problem expressed in the form of equation (2) is to be solved, equation (2) is often written in the incremental form, and the z transform [17] can be used to explore an algorithm’s stability for the cases of stiffness hardening and softening. The solution to equation (2) is also often written in the incremental form, as in the following equation:(50)mat+Δt−at+cvt+Δt−vt+ktanut+Δt−ut=ft+Δt−ft,where ktan is the tangent stiffness at time step t+Δt [18].When addressing structural dynamic nonlinearity in the form of equation (2), the CQ-2x method uses equations (46), (21), and (50).The stability of an integration algorithm is generally analyzed using one of two approaches, namely, recurrence relations using an amplification matrix or discrete control theory. The unconditional stability of CQ-2x when solving nonlinear equation (3) has been analyzed through the use of an amplification matrix. In this section, discrete control theory will be used to analyze the stability of CQ-2x when solving nonlinear equation (2). The discrete z transform and its real translation property are defined as follows [18]:(51)Xz=Zxn=∑n=0∞xnz−n,Zxn−k=z−kXz,Zui+1Δt=Xz,Zvi+1Δt=Xvz,Zai+1Δt=Xaz,Zfi+1Δt=Fz.The stability of CQ-2x when solving nonlinear equation (2) can be analyzed by applying discrete control theory to equations (46), (21), and (50).First, the stability analysis of equation (46) using the z transform is considered. By applying the z transform and its real translation property to equation (46), this equation can be solved to determine X(z) in terms of the system properties, the time step, and F(z). This solution leads to the discrete transfer function of the algorithm, which is defined as the ratio of the z transforms of the output ui+1 and the input fi+1 and can be expressed in the following general form:(52)XzFz=4zd2z2+d1z+d0,where(53)d2=mx2+2x+1Ωt+Δt2+4ξt+ΔtΩt+Δt+4,d1=m2−2x2Ωt+Δt2−8,d0=mx2−2x+1Ωt+Δt2−4ξt+ΔtΩt+Δt+4.It can be seen that the denominator of the transfer function in equation (52) has the same form as the characteristic equation given in equation (C1) for the amplification matrix of equation (46), whose unconditional stability has been proven.Since the displacement calculation is unconditionally stable, the velocity calculation in equation (21), which depends solely on the displacement, is also naturally unconditionally stable.Now, the stability of equation (50) will be analyzed. By combining equations (46), (21), and (50), the transfer function defined as the ratio of the z transforms of the output ai+1 and the input fi+1 can be expressed in the following general form:(54)XazFz=n3z3+n2z2+n1z+n0d3z3+d2z2+d1z,where(55)d3=mΔtx2+2x+1Ωt+Δt2+4ξt+ΔtΩt+Δt+4,d2=mΔt2−2x2Ωt+Δt2−8,d1=mΔtx2−2x+1Ωt+Δt2−4ξt+ΔtΩt+Δt+4.Sincen3, n2, n1, and n0 do not affect the algorithm’s stability, for brevity, they are not listed in this paper.Although the denominators of equations (52) and (54) differ by a multiplicative factor Δt, their representative characteristic equations are the same, which means that the acceleration calculation using equation (50) is also unconditionally stable. An analysis shows that the ratio of ktan to kt+Δt only affects the value of the numerator of equation (54). Thus, it can be concluded that CQ-2x is unconditionally stable when it is used to solve the structural nonlinear dynamics equation given in equation (2); this will be validated by Example 2.


### 5.5. Accuracy Analysis for the Linear System

The accuracy of an integration algorithm generally depends on its numerical dispersion and energy dissipation characteristics. For the purposes of an accuracy analysis, the solutions to a free vibration problem obtained using CQ-2x can be written in the form of equation (D3). The two principal roots (no spurious real root) can be used to analyze the numerical dispersion and energy dissipation characteristics. The first characteristic is generally expressed in terms of the relative period error, as defined in equation (D6). The latter one can be expressed in terms of the numerical damping ratio ξ¯, as defined in equation (D4).Figures5 and 6 show the relative period errors (Pe) and the equivalent damping ratios ξ¯, respectively, of the proposed CQ-2x method and the KR-α method for various values of ρ∞. As shown in the figures, both Pe and ξ¯ increase with decreasing ρ∞, and ρ∞=0 produces the maximum Pe and ξ¯.Figure 5
Relative period error of CQ-2x when ξ=0.Figure 6
Numerical damping ratio of CQ-2x when ξ=0.Figure5 shows that the numerical dispersion of CQ-2x is similar to that of the KR-α method for the same value of ρ∞ when ρ∞≠1. Figure 6 shows that the KR-α method’s energy dissipation is better than that of CQ-2x. However, at small values of Ω, both Pe and ξ¯ are small for CQ-2xρ∞∈0.8,1 regardless of the value of ρ∞, which means that the lower-mode responses in an MDOF system are negligibly affected by the numerical damping. By contrast, ξ¯ increases with increasing Ω, indicating that the undesired higher modes can be effectively damped out. Note that, for ρ∞=1, the CQ-2x method exhibits no numerical energy dissipation, and the Pe value becomes identical to that of the KR-αρ∞=1 method and Newmark’s CAA method.


### 5.6. Self-Starting and Overshoot

There are several possible methods of self-starting, which is a term that refers to a parameterless unified starting procedure [26]. In this paper, the self-starting of CQ-2x relies on u−Δt, as calculated using the CDM, as shown in the following equation:(56)u−Δt=u0−vΔt+12a0Δt2.Whent=−Δt, equation (46) can be written as(57)m+12cΔt+x+124k0Δt2uΔt=2m+x2−12k0Δt2u0+−m+12cΔt−x−124k0Δt2u−Δt+f0Δt2.Substituting equation (56) into equation (57) yields a value for uΔt as follows:(58)m+12cΔt+x+124k0Δt2uΔt=m+12cΔt+x2+2x−34k0Δt2u0+m−12cΔt+x−124k0Δt2Δtv0+12−m+12cΔt−x−124k0Δt2Δt2a0+f0Δt2.With known initial values ofu0, v0, and a0, the value of uΔt can be calculated using equation (58). Thus, a self-starting displacement calculation procedure for CQ-2x can be obtained.When the velocity is to be calculated, equation (59) can be used, with a second-order time accuracy of velocity:(59)vΔt=3uΔt−3u0Δt−2v0−12a0Δt.Equation (59) represents the self-starting velocity calculation in CQ-2x.WhenΩ0⟶∞, equation (58) becomes(60)limΩ0⟶∞uΔt=x2+2x−3x+12u0+x−12x+12Δtv0−12x−12x+12Δt2a0.For CQ-2x, as expressed in equation (46), the recursive relationship among ut+2Δt, ut+Δt, and ut can be expressed as follows:(61)limΩ⟶∞ut+2Δt=2x−11+xut+Δt−1−x21+x2ut.It can be seen from equations (60) and (61) that CQ-2x has zero-order displacement overshoot. Similarly, it is easy to find that CQ-2x has zero-order velocity overshoot by considering equations (21) and (59). Thus, CQ-2x can be classified as a family of noniterative algorithms of the form U0V0ρ∞,ρ∞ [20].


### 5.7. Calculation Procedure for MDOF Systems

When applied in MDOF systems, the forms of the aforementioned equations have some changes. When applied in MDOF systems,m,ct+iΔt,andkt+iΔt appearing in a SDOF system should then be replaced by matrices M,Ct+iΔt,andKt+iΔt; ut+iΔt,vt+iΔt,at+iΔt,andft+iΔt appearing in a SDOF system should be replaced by vectors ut+iΔt,vt+iΔt,at+iΔt,andft+iΔt. The calculation procedure for MDOF systems is provided as follows:(1)
Formulate the stiffness matrixK0, mass matrix M, and damping matrix C0 of the systems.(2)
Select appropriate coefficientρ∞ and time increment Δt. Determine x by solving equation (49d).(3)
Obtain initial valuesu0andv0. Calculate a0 by(62)a0=M−1f0−C0v0−K0u0.(4)
CalculateuΔtandvΔt by equations (58) and (59), respectively. Then, calculate nonlinear stiffness and damping KΔtandCΔt by equation (4). Subsequently, calculate aΔt by the following equation:(63)at+Δt=1mft+Δt−ct+Δtvt+Δt−kt+Δtut+Δt.(5)
For each time step (i = 2, 4, 5, …, n), calculate uiΔt and viΔt by equations (49a) and (49b), respectively. Then, calculate nonlinear stiffness and damping KiΔtandCiΔt by equation (4). Subsequently, calculate aiΔt by equation (49c).The main calculation efforts of CQ-2x are concentrated in Step 5 discussed above. The main calculation efforts in Step 5 are concentrated in equation (49a) because the coefficient matrix (M+1/2Ct+ΔtΔt+x+12/4Kt+ΔtΔt2) in equation (49a) is nondiagonal and solving the coupling linear equations needs ϑN3 times of multiplication or division. The computational costs consumed by equation (49b) are ϑN, and those consumed by equations (4) and (3) are ϑN2 since solving equation (4) only needs vector calculations and solving equation (3) only needs solving linear decoupling equations whose coefficient matrix is the diagonal matrix M. Thus, it is true that there is only one solver of equation (49a) within one time step for CQ-2x. Furthermore, when no velocity and acceleration need to be calculated, the entire calculation of CQ-2x can be completed with only the involvement of displacement. In this case, the time consumption of CQ-2x will become less, and the storage requirements of CQ-2x will become considerably lower than other unconditionally stable noniterative algorithms.


## 5.1. Determination ofb7

The unconditional stability of an NDA of the form given in equation (20) means that when Ω∈0,+∞ and ξ∈0,+∞, equations (32a)–(32c) are always satisfied, which requires that the following must hold:(34)1−b7≥0,1+b8+b5≥0,2−4b7≥0,1+2b5≥0,1+2b8+b5≥0.An intermediate variablex is introduced:(35)x=1+2b8+b5.Consequently,b8 can be written as(36)b8=12x−b5−1.Then, simplifying equation (34) leads to(37)b7≤0.5,b8≤x,b5≥−0.5,x≥0.An analysis shows that, for a damping structure, a largerb7 leads to more satisfactory algorithmic characteristics; thus, b7 is assigned its largest possible value of 0.5 according to equation (37):(38)b7=0.5.Substituting equations (36) and (38) into equation (20) yields the preliminary algorithmic form of the CQ-2x method:(39)m+0.5ct+ΔtΔt+12x+b5+1kt+ΔtΔt2ut+2Δt=2m+b5kt+ΔtΔt2ut+Δt+−m+0.5ct+ΔtΔt+12x−b5−1kt+ΔtΔt2ut+ft+ΔtΔt2.It can be seen from the above analysis that there are only two coefficients,b5 and x, that are undetermined in the preliminary form of the CQ-2x method. These coefficients will be determined by introducing controllable numerical dissipation into CQ-2x.


## 5.2. Determination ofb5 and x

Whenξ=0, according to equation (36), equations (31a) and (31b) can be written as(40a)%A1=A11+A222=1+1/2b5Ω21+1/2b5+1/2x+1/2Ω2,(40b)A2=1+1/2+1/2b5−1/2xΩ21+1/2+1/2b5+1/2xΩ2.The effects ofb5 and x on the relative period error of the preliminary CQ-2x method have been analyzed. The results when Δt/T=0.1, b5=−0.5,−0.3,0,0.3,1,2,and4, and x∈0,4 are plotted in Figure 2. It can be seen that, as b5 and x increase, the relative period error also increases. Thus, to achieve the best accuracy, x and b5 should be assigned their lowest possible values of 0 and −0.5, respectively. The analysis shows that when x=0 and b5=−0.5, the preliminary CQ-2x method has a nominal third-order time accuracy of displacement and a global second-order time accuracy; this will be validated by Example 1. However, if numerical dissipation is needed to filter out spurious vibrations in high-vibration modes, neither x=0 nor b5=−0.5 is appropriate.Figure 1
NDA stability domain.Figure 2
Variations in the relative period error of the preliminary CQ-2x method with b5 and x when Δt/T=0.1.It can be proven that whenx=0, the numerical damping ratio is also ξ¯=0, as shown in the following equation:(41)x=0,A2≥A12,A2=1+1/2+1/2b5Ω21+1/2+1/2b5Ω2=1,ρA=A2=1,ξ¯=0.It can also be proven that whenb5=−0.5, the numerical damping ratio ξ¯ of very high vibration modes Ω⟶∞ is also zero, as shown in the following equation:(42)b5=−0.5,A2≤A12,limΩ⟶∞ρA=limΩ⟶∞A1−A12−A2=−1/4−1/2x1/2x+1/4=1,limΩ⟶∞ξ¯=0.To introduce an appropriate level of numerical dissipation,b5 and x must be greater than their lowest values of 0 and −0.5, respectively. The numerical damping ratio depends on the spectral radius of the amplification matrix. It may be preferable for the spectral radius of the amplification matrix A of the preliminary CQ-2x method to decrease with increasing Ω with no bifurcation. According to equations (D1a) and (D1b), no bifurcation occurs when A2≥A12, which requires the following:(43)A2−A12=−x24−b52−14Ω4+Ω2≥0,whenΩ∈0,+∞.Equation (43) is satisfied when(44)b5≥x22−12.Since a lower value ofb5 results in a higher accuracy, as mentioned before, b5 is assigned the lowest possible value according to equation (44):(45)b5=x22−12.Substituting equation (45) into equation (39) yields(46)m+12ct+ΔtΔt+x+124kt+ΔtΔt2ut+2Δt=2m+x2−12kt+ΔtΔt2ut+Δt+−m+12ct+ΔtΔt−x−124kt+ΔtΔt2ut+ft+ΔtΔt2.Equation (46) is the final algorithmic form of displacement calculation in the CQ-2x method.To determinex, while simultaneously confirming the unconditional stability of the CQ-2x method, consider the following two extreme cases: Ω⟶0 and Ω⟶∞. In accordance with equations (31a) and (31b), (36), (45), and (46) and considering the limits, it can be shown that equation (C1) reduces to the following:(47a)λ−12=0,whenΩ⟶0,(47b)λ−x−1x+12=0,whenΩ⟶∞.It can be calculated from equation (47b) that(48)x=1−ρ∞1+ρ∞.Figures3 and 4 show how the spectral radius of the CQ-2x method varies with ρ∞ and Ω; the findings are fully consistent with the above analysis.Figure 3
Spectral radius of CQ-2x when ξ=0.Figure 4
Spectral radius of CQ-2x when ξ=0.05.


## 5.3. The Proposed Algorithm

Now that all coefficients in the preliminary NDA form presented in equation (6) have been determined, a family of unconditionally stable noniterative algorithms, collectively referred to as the CQ-2x method with one-parameter x, can be obtained by combining equations (21), (23), (46), and (48), which are rewritten as(49a)m+12ct+ΔtΔt+x+124kt+ΔtΔt2ut+2Δt=2m+x2−12kt+ΔtΔt2ut+Δt+−m+12ct+ΔtΔt−x−124kt+ΔtΔt2ut+ft+ΔtΔt2,(49b)vt+2Δt=3ut+2Δt−4ut+Δt+ut2Δt,(49c)at+2Δt=1mft+2Δt−ct+2Δtvt+2Δt−kt+2Δtut+2Δt,where(49d)x=1−ρ∞1+ρ∞.


## 5.4. Stability for Systems with Nonlinear Stiffness

It has been demonstrated that the CQ-2x method is unconditionally stable when it is used to solve the nonlinear structural dynamics equation given in equation (3). Here, the stability of CQ-2x when it is used to solve the nonlinear structural dynamics equation given in equation (2) will be analyzed. When a general nonlinear structural dynamics problem expressed in the form of equation (2) is to be solved, equation (2) is often written in the incremental form, and the z transform [17] can be used to explore an algorithm’s stability for the cases of stiffness hardening and softening. The solution to equation (2) is also often written in the incremental form, as in the following equation:(50)mat+Δt−at+cvt+Δt−vt+ktanut+Δt−ut=ft+Δt−ft,where ktan is the tangent stiffness at time step t+Δt [18].When addressing structural dynamic nonlinearity in the form of equation (2), the CQ-2x method uses equations (46), (21), and (50).The stability of an integration algorithm is generally analyzed using one of two approaches, namely, recurrence relations using an amplification matrix or discrete control theory. The unconditional stability of CQ-2x when solving nonlinear equation (3) has been analyzed through the use of an amplification matrix. In this section, discrete control theory will be used to analyze the stability of CQ-2x when solving nonlinear equation (2). The discrete z transform and its real translation property are defined as follows [18]:(51)Xz=Zxn=∑n=0∞xnz−n,Zxn−k=z−kXz,Zui+1Δt=Xz,Zvi+1Δt=Xvz,Zai+1Δt=Xaz,Zfi+1Δt=Fz.The stability of CQ-2x when solving nonlinear equation (2) can be analyzed by applying discrete control theory to equations (46), (21), and (50).First, the stability analysis of equation (46) using the z transform is considered. By applying the z transform and its real translation property to equation (46), this equation can be solved to determine X(z) in terms of the system properties, the time step, and F(z). This solution leads to the discrete transfer function of the algorithm, which is defined as the ratio of the z transforms of the output ui+1 and the input fi+1 and can be expressed in the following general form:(52)XzFz=4zd2z2+d1z+d0,where(53)d2=mx2+2x+1Ωt+Δt2+4ξt+ΔtΩt+Δt+4,d1=m2−2x2Ωt+Δt2−8,d0=mx2−2x+1Ωt+Δt2−4ξt+ΔtΩt+Δt+4.It can be seen that the denominator of the transfer function in equation (52) has the same form as the characteristic equation given in equation (C1) for the amplification matrix of equation (46), whose unconditional stability has been proven.Since the displacement calculation is unconditionally stable, the velocity calculation in equation (21), which depends solely on the displacement, is also naturally unconditionally stable.Now, the stability of equation (50) will be analyzed. By combining equations (46), (21), and (50), the transfer function defined as the ratio of the z transforms of the output ai+1 and the input fi+1 can be expressed in the following general form:(54)XazFz=n3z3+n2z2+n1z+n0d3z3+d2z2+d1z,where(55)d3=mΔtx2+2x+1Ωt+Δt2+4ξt+ΔtΩt+Δt+4,d2=mΔt2−2x2Ωt+Δt2−8,d1=mΔtx2−2x+1Ωt+Δt2−4ξt+ΔtΩt+Δt+4.Sincen3, n2, n1, and n0 do not affect the algorithm’s stability, for brevity, they are not listed in this paper.Although the denominators of equations (52) and (54) differ by a multiplicative factor Δt, their representative characteristic equations are the same, which means that the acceleration calculation using equation (50) is also unconditionally stable. An analysis shows that the ratio of ktan to kt+Δt only affects the value of the numerator of equation (54). Thus, it can be concluded that CQ-2x is unconditionally stable when it is used to solve the structural nonlinear dynamics equation given in equation (2); this will be validated by Example 2.


## 5.5. Accuracy Analysis for the Linear System

The accuracy of an integration algorithm generally depends on its numerical dispersion and energy dissipation characteristics. For the purposes of an accuracy analysis, the solutions to a free vibration problem obtained using CQ-2x can be written in the form of equation (D3). The two principal roots (no spurious real root) can be used to analyze the numerical dispersion and energy dissipation characteristics. The first characteristic is generally expressed in terms of the relative period error, as defined in equation (D6). The latter one can be expressed in terms of the numerical damping ratio ξ¯, as defined in equation (D4).Figures5 and 6 show the relative period errors (Pe) and the equivalent damping ratios ξ¯, respectively, of the proposed CQ-2x method and the KR-α method for various values of ρ∞. As shown in the figures, both Pe and ξ¯ increase with decreasing ρ∞, and ρ∞=0 produces the maximum Pe and ξ¯.Figure 5
Relative period error of CQ-2x when ξ=0.Figure 6
Numerical damping ratio of CQ-2x when ξ=0.Figure5 shows that the numerical dispersion of CQ-2x is similar to that of the KR-α method for the same value of ρ∞ when ρ∞≠1. Figure 6 shows that the KR-α method’s energy dissipation is better than that of CQ-2x. However, at small values of Ω, both Pe and ξ¯ are small for CQ-2xρ∞∈0.8,1 regardless of the value of ρ∞, which means that the lower-mode responses in an MDOF system are negligibly affected by the numerical damping. By contrast, ξ¯ increases with increasing Ω, indicating that the undesired higher modes can be effectively damped out. Note that, for ρ∞=1, the CQ-2x method exhibits no numerical energy dissipation, and the Pe value becomes identical to that of the KR-αρ∞=1 method and Newmark’s CAA method.


## 5.6. Self-Starting and Overshoot

There are several possible methods of self-starting, which is a term that refers to a parameterless unified starting procedure [26]. In this paper, the self-starting of CQ-2x relies on u−Δt, as calculated using the CDM, as shown in the following equation:(56)u−Δt=u0−vΔt+12a0Δt2.Whent=−Δt, equation (46) can be written as(57)m+12cΔt+x+124k0Δt2uΔt=2m+x2−12k0Δt2u0+−m+12cΔt−x−124k0Δt2u−Δt+f0Δt2.Substituting equation (56) into equation (57) yields a value for uΔt as follows:(58)m+12cΔt+x+124k0Δt2uΔt=m+12cΔt+x2+2x−34k0Δt2u0+m−12cΔt+x−124k0Δt2Δtv0+12−m+12cΔt−x−124k0Δt2Δt2a0+f0Δt2.With known initial values ofu0, v0, and a0, the value of uΔt can be calculated using equation (58). Thus, a self-starting displacement calculation procedure for CQ-2x can be obtained.When the velocity is to be calculated, equation (59) can be used, with a second-order time accuracy of velocity:(59)vΔt=3uΔt−3u0Δt−2v0−12a0Δt.Equation (59) represents the self-starting velocity calculation in CQ-2x.WhenΩ0⟶∞, equation (58) becomes(60)limΩ0⟶∞uΔt=x2+2x−3x+12u0+x−12x+12Δtv0−12x−12x+12Δt2a0.For CQ-2x, as expressed in equation (46), the recursive relationship among ut+2Δt, ut+Δt, and ut can be expressed as follows:(61)limΩ⟶∞ut+2Δt=2x−11+xut+Δt−1−x21+x2ut.It can be seen from equations (60) and (61) that CQ-2x has zero-order displacement overshoot. Similarly, it is easy to find that CQ-2x has zero-order velocity overshoot by considering equations (21) and (59). Thus, CQ-2x can be classified as a family of noniterative algorithms of the form U0V0ρ∞,ρ∞ [20].


## 5.7. Calculation Procedure for MDOF Systems

When applied in MDOF systems, the forms of the aforementioned equations have some changes. When applied in MDOF systems,m,ct+iΔt,andkt+iΔt appearing in a SDOF system should then be replaced by matrices M,Ct+iΔt,andKt+iΔt; ut+iΔt,vt+iΔt,at+iΔt,andft+iΔt appearing in a SDOF system should be replaced by vectors ut+iΔt,vt+iΔt,at+iΔt,andft+iΔt. The calculation procedure for MDOF systems is provided as follows:(1)
Formulate the stiffness matrixK0, mass matrix M, and damping matrix C0 of the systems.(2)
Select appropriate coefficientρ∞ and time increment Δt. Determine x by solving equation (49d).(3)
Obtain initial valuesu0andv0. Calculate a0 by(62)a0=M−1f0−C0v0−K0u0.(4)
CalculateuΔtandvΔt by equations (58) and (59), respectively. Then, calculate nonlinear stiffness and damping KΔtandCΔt by equation (4). Subsequently, calculate aΔt by the following equation:(63)at+Δt=1mft+Δt−ct+Δtvt+Δt−kt+Δtut+Δt.(5)
For each time step (i = 2, 4, 5, …, n), calculate uiΔt and viΔt by equations (49a) and (49b), respectively. Then, calculate nonlinear stiffness and damping KiΔtandCiΔt by equation (4). Subsequently, calculate aiΔt by equation (49c).The main calculation efforts of CQ-2x are concentrated in Step 5 discussed above. The main calculation efforts in Step 5 are concentrated in equation (49a) because the coefficient matrix (M+1/2Ct+ΔtΔt+x+12/4Kt+ΔtΔt2) in equation (49a) is nondiagonal and solving the coupling linear equations needs ϑN3 times of multiplication or division. The computational costs consumed by equation (49b) are ϑN, and those consumed by equations (4) and (3) are ϑN2 since solving equation (4) only needs vector calculations and solving equation (3) only needs solving linear decoupling equations whose coefficient matrix is the diagonal matrix M. Thus, it is true that there is only one solver of equation (49a) within one time step for CQ-2x. Furthermore, when no velocity and acceleration need to be calculated, the entire calculation of CQ-2x can be completed with only the involvement of displacement. In this case, the time consumption of CQ-2x will become less, and the storage requirements of CQ-2x will become considerably lower than other unconditionally stable noniterative algorithms.


## 6. Numerical Examples

To evaluate the overall behavior and confirm the convergence rate, controllable numerical dissipation, unconditional stability in case of stiffness hardening, and computational efficiency of the derived algorithms, four examples are given. The CAA method, as an example of an implicit algorithm, and the KR-α method, as an example of a family of unconditionally stable structure-dependent noniterative algorithms, are considered for comparison.


### 6.1. Example 1

The purpose of this example is to validate the convergence rate of the CQ-2x method. Consider a single-DOF linear resonance problem with the initial conditions u0=v0=1:(64)at+ω2ut=sinωt.The exact solutions [27] to this problem are(65a)ut=2ω+12ω2sinωt+2ω−t2ωcosωt,(65b)vt=t−2ω2sinωt+cosωt,(65c)at=1−2ω2sinωt+ωt−2ω2cosωt.The circular self-frequencyω was assigned a value of 2π. Thus, the period is T=1 s. The integration interval was chosen to be TFin=10T=10 s. In the calculations, the step size Δt was varied from 10−1 s to 10−4 s. The plots in Figures 7–12 show the curves of the maximum relative absolute errors and average relative absolute errors in displacement, velocity, and acceleration for the duration of the simulation. These errors were computed as functions of Δt (on a logarithmic scale) using the CQ-2x, CAA, and KR-α methods. When ρ∞≠1, the CQ-2x method shows first-order rates of convergence in displacement, velocity, and acceleration, identical to the KR-α method. When ρ∞=1, CQ-2x shows second-order rates of convergence in displacement, velocity, and acceleration, identical to the CAA method.Figure 7
Convergence rate of CQ-2x assessed in terms of the maximum relative displacement error.Figure 8
Convergence rate of CQ-2x assessed in terms of the average relative displacement error.Figure 9
Convergence rate of CQ-2x assessed in terms of the maximum relative velocity error.Figure 10
Convergence rate of CQ-2x assessed in terms of the average relative velocity error.Figure 11
Convergence rate of CQ-2x assessed in terms of the maximum relative acceleration error.Figure 12
Convergence rate of CQ-2x assessed in terms of the average relative acceleration error.


### 6.2. Example 2

The purpose of this example is to validate the controllable numerical dissipation and unconditional stability in case of stiffness hardening of the CQ-2x method. In this example, the well-known unforced and undamped Duffing oscillator [28] is considered. The nonlinear dynamics equation is(66)a+S1u1+S2u2=0,where S1=100 and S2−10. The initial conditions are u0=1.5 and v0=0.Equation (66) represents a hardening elastic spring. This is an example of a conservative system that maintains a constant total energy E in its exact solutions. For equation (66),(67)E=12v2+12S1u2+14S1S2u4.To assess the time integration scheme, the percentage error in terms of the energy is introduced [29]:(68)Er=E−E0E0×100%,where E0 is the total energy at t=0.Equation (66) has an exact solution expressed in terms of Jacobian elliptic functions [30]:(69)ut=u0cn−w^t,k,where(70)ω^2=S11+S2u02,k2=S2u0221+S2u02.The periodT of this oscillator is equal to 0.15. The exact phase portrait is given in Figure 13. Also, shown in Figures 14–22 are the phase portraits of numerical results obtained using different integration schemes with the same time step size of Δt=T/20 over a time duration of 100 T. The maximum values of the percentage errors of the nine considered algorithms over this time duration of 100 T are summarized in Table 1 for various time step sizes. When Δt=1/25T, CQ-2xρ∞=1 has an error of 6.06%, whereas the KR-α method has an error of at least 99.63%. It can be seen from Figures 14–17 that ρ∞ is smaller, the numerical damping is bigger, and the undamped Duffing oscillator decays faster. To validate the unconditional stability of CQ-2x in the stiffness hardening case, a time step size of Δt=1/2T over a time duration of 100 T is considered. Figures 23–31 show the phase portraits of the numerical results obtained using the nine integration schemes with the same time step size of Δt=1/2T over a time duration of 100 T. It is meaningless to study an algorithm’s accuracy when Δt=1/2T, but it is useful to explore its stability. Figures 23–26 show that even when Δt=1/2T, CQ-2xρ∞=1,0.8,0.5,0 is still stable. By contrast, the KR-α method is robustly stable only when ρ∞=0, as seen in Figure 30; when ρ∞=1,0.8,and0.5, the KR-α method becomes evidently unstable, as seen in Figures 27–29. e CAA method also exhibits some instability, as seen in Figure 31 (for the CAA method, the iteration termination criterion is that the absolute value of the error between the displacements in the previous and current iterations is less than 10−6).Figure 13
Exact solution for the phase portrait of the hardening spring.Figure 14
Phase portrait of the hardening spring calculated using CQ-2xρ∞=1 when Δt=T/20.Figure 15
Phase portrait of the hardening spring calculated using CQ-2xρ∞=0.8 when Δt=T/20.Figure 16
Phase portrait of the hardening spring calculated using CQ-2xρ∞=0.5 when Δt=T/20.Figure 17
Phase portrait of the hardening spring calculated using CQ-2xρ∞=0 when Δt=T/20.Figure 18
Phase portrait of the hardening spring calculated using the KR-α method ρ∞=1 when Δt=T/20.Figure 19
Phase portrait of the hardening spring calculated using the KR-α method ρ∞=0.8 when Δt=T/20.Figure 20
Phase portrait of the hardening spring calculated using the KR-α method ρ∞=0.5 when Δt=T/20.Figure 21
Phase portrait of the hardening spring calculated using the KR-α method ρ∞=0 when Δt=T/20.Figure 22
Phase portrait of the hardening spring calculated using the CAA method whenΔt=T/20.Table 1
Maximum percentage errors of the total energy over a time duration of 100T for the hardening spring.
CQ-2xKR-αCAAΔtρ∞=1ρ∞=0.8ρ∞=0.5ρ∞=0ρ∞=1ρ∞=0.8ρ∞=0.5ρ∞=0T/10000.0019.8845.5477.2271.6368.5560.601.310.00T/2000.0960.8087.9398.2496.0095.2292.976.830.04T/1000.3679.7696.0799.6899.3298.2597.3615.890.17T/501.4792.1399.0899.97158.3999.3599.0437.140.66T/256.0697.7399.87100.00105670.0299.7499.63400.592.62T/209.5198.6099.94100.00∞99.8199.73700089.874.07T/1516.7799.2899.98100.00∞99.8799.81219047.467.13T/2FiniteFiniteFiniteFinite∞∞∞finite∞Figure 23
Phase portrait of the hardening spring calculated using CQ-2xρ∞=1 when Δt=1/2T.Figure 24
Phase portrait of the hardening spring calculated using CQ-2xρ∞=0.8 when Δt=1/2T.Figure 25
Phase portrait of the hardening spring calculated using CQ-2xρ∞=0.5 when Δt=1/2T.Figure 26
Phase portrait of the hardening spring calculated using CQ-2xρ∞=0 when Δt=1/2T.Figure 27
Phase portrait of the hardening spring calculated using the KR-α method ρ∞=1 when Δt=1/2T.Figure 28
Phase portrait of the hardening spring calculated using the KR-methodρ∞=0.8 when Δt=1/2T.Figure 29
Phase portrait of the hardening spring calculated using the KR-α method ρ∞=0.5 when Δt=1/2T.Figure 30
Phase portrait of the hardening spring calculated using the KR-α method ρ∞=0 when Δt=1/2T.Figure 31
Phase portrait of the hardening spring calculated using the CAA method whenΔt=1/2T.The time consumptions of the different algorithms with different time steps over the time duration of 100T are listed in Table 2. For the same time step size, the noniterative algorithms (CQ-2x and KR-α) are more efficient than the implicit algorithm (CAA). In addition, the CQ-2x method is more efficient than the KR-α method for the same time step size. For example, when Δt=1/25T, the time consumption of CQ-2x is at most 0.0015 seconds, whereas the KR-α method requires at least 0.0048 seconds.Table 2
Time consumptions of the different algorithms for solving the Duffing vibration problem in Example 2 over a time duration of 100T (unit: second).
CQ-2xKR-αCAAΔtρ∞=1ρ∞=0.8ρ∞=0.5ρ∞=0ρ∞=1ρ∞=0.8ρ∞=0.5ρ∞=0T/10000.05900.00900.01800.01740.33050.18710.19360.19471.1631T/2000.01150.00180.00370.00360.06860.04050.04740.04310.2700T/1000.00600.00090.00200.00190.03310.01990.01980.02410.1386T/500.00320.00040.00090.00090.01660.00950.01040.01250.0835T/250.00150.00020.00050.00050.00790.00480.00480.00800.0511T/200.00150.00020.00040.00040.00660.00400.00390.00680.0426T/150.00090.00010.00030.00030.00580.00270.00290.00510.0386


### 6.3. Example 3

The purpose of this example is to validate the efficiency of the CQ-2x method compared with the implicit method and KR−α method. In this example, the forced vibration response of a n-DOF (n = 200, 500, 1000, 2000, and 4000) spring-mass system was considered, as shown in Figure 32. The structural properties of this system are assumed to be mi=100kg and ki=1071−ui−ui−12N/m, in which i=1,2,…,n. This system is subjected to a ground acceleration of 10sinπt. The highest natural frequency is 632.4 rad/s [31]. According to the highest natural frequency, a time step duration of Δt=0.001s is obtained from Newmark’s average acceleration method and is selected to be the exact solution. The responses of displacement and velocity versus time of the 200th DOF when N = 200 are plotted in Figures 33 and 34 . The analysis is accomplished by considering the time step duration equal to 0.02 s. The maximum displacement and velocity error ratios of the 200th DOF when n = 200 are listed in Table 3. It can be seen that the accuracy of the CQ-2x method is sufficient. The maximum error ratio here is defined by the function of max (abs(calculated value-exact value)))/max(abs(exact value)), where max = maximum and abs = absolute. Importantly, without the need for iteration calculations, CQ-2x algorithms have approximately less than half the computational time consumption compared to that consumed by the implicit algorithm CAA. With solving the inverse matrix only once for each time step, CQ-2x algorithms are more efficient than the noniterative KR−α method in which solving the inverse matrix twice for each time step is needed, as listed in Table 4.Figure 32
Schematic plot of then-DOF spring-mass system.Figure 33
Displacement responses of the 200th DOF whenn  = 200.Figure 34
Velocity responses of the 200th DOF whenn  = 200.Table 3
Maximum displacement and velocity responses’ errors of the 200th DOF whenn = 200.
CAACQ-2xKR-αρ∞=1ρ∞=0.8ρ∞=1ρ∞=0.8Displacement0.003310.002840.031170.002830.00289Velocity0.006030.005530.034540.029550.03537Table 4
Time consumptions of the different algorithms for solving then-DOF system in Example 3 (unit: second).
CAACQ-2xKR-αρ∞=1ρ∞=0.8ρ∞=1ρ∞=0.8n=2005.374781.633621.516273.752333.54338n=50039.1792314.8079711.1223131.6716823.96452n=1000246.3178986.0823991.89233182.75176183.34749n=20002045.64400318.19379492.85635893.826941094.05923n=400014992.5984811.7546220.65112205.04014868.966


### 6.4. Example 4

The purpose of this example is to validate the efficiency of the CQ-2x method in case of Vehicle-Bridge Coupled Vibration where the contact nonlinearity exists. In this example, a simple beam of span length = 25 m is subjected to a moving spring mass [32] (as shown in Figure 35). The following data are adopted: Young’s modulus E = 2.87 Gpa, ρ=0.2, moment of inertia I=2.90m4, mass per unit length m=2303kg/m, suspended mass Mv=5750kg, suspension stiffness kv=1595kN/m, and speed v=90km/h. The frequency of vibration computed of the spring mass is ωv=16.66rad/s. The frequencies of vibration computed of the bridge are ω1=30.02rad/s,ω2=120.08rad/s,ω3=270.19rad/s,ω4=480.34rad/s,ω5=750.58rad/s,andω6=1081.0rad/s. The maximum frequency of vibration computed of the bridge is ωmax=95432rad/s. According to stability requirement of the CDM, the widest time step satisfying CDM’s stability is 2/ωmax=2.0957×10−5s. A time step duration of Δt=2×10−5s is obtained for the CDM, and the calculation result is selected to be the exact solution, while with unconditional stability, the time step duration of CQ2xρ∞=1, KR-αρ∞=1, and CAA is only determined by the requirements of the calculation accuracy; here, Δt=0.001s. The midpoint vertical deflection of the beam and deflection of the spring mass versus time are plotted in Figures 36 and 37. The maximum displacement error ratios (defined as that in Example 3) of the midpoint vertical deflection of the beam and the deflection of the spring mass are listed in Table 5. It can be seen that the accuracy of the CQ-2x method is sufficient. The efficiency of the CDM Δt=0.00002s, CQ2xρ∞=1, KR-αρ∞=1, and CAA Δt=0.001s is compared. The average time consumption of the four algorithms in which every algorithm had been performed five times is listed in Table 6. In this example, limited by its conditional stability, the CDM has to be performed with very small Δt=0.00002s, which leads to its efficiency much lower than the other three algorithms. Without the need for iteration calculations, CQ-2x and KR-α algorithms show higher efficiency than that of the implicit algorithm CAA although the maximum number of iteration in each time step of CAA is only 3. The CQ-2x algorithm are more efficient than the noniterative KR−α algorithm since the computation amount in each time step of CQ-2x is much less than that of KR-α.Figure 35
Beam with the moving spring mass.Figure 36
Deflection of the spring mass.Figure 37
Midpoint vertical deflection of the beam.Table 5
Maximum deflection error ratios in Example 4.
CAACQ-2xρ∞=1KR-αρ∞=1Midpoint of beam0.06720.01770.0061Spring mass0.03540.03030.0051Table 6
Time consumption of the different algorithms for solving the beam with the moving spring mass in Example 4 (unit: second).
CDMCAACQ-2xρ∞=1KR-αρ∞=11153.891.409380.512521.21876


## 6.1. Example 1

The purpose of this example is to validate the convergence rate of the CQ-2x method. Consider a single-DOF linear resonance problem with the initial conditions u0=v0=1:(64)at+ω2ut=sinωt.The exact solutions [27] to this problem are(65a)ut=2ω+12ω2sinωt+2ω−t2ωcosωt,(65b)vt=t−2ω2sinωt+cosωt,(65c)at=1−2ω2sinωt+ωt−2ω2cosωt.The circular self-frequencyω was assigned a value of 2π. Thus, the period is T=1 s. The integration interval was chosen to be TFin=10T=10 s. In the calculations, the step size Δt was varied from 10−1 s to 10−4 s. The plots in Figures 7–12 show the curves of the maximum relative absolute errors and average relative absolute errors in displacement, velocity, and acceleration for the duration of the simulation. These errors were computed as functions of Δt (on a logarithmic scale) using the CQ-2x, CAA, and KR-α methods. When ρ∞≠1, the CQ-2x method shows first-order rates of convergence in displacement, velocity, and acceleration, identical to the KR-α method. When ρ∞=1, CQ-2x shows second-order rates of convergence in displacement, velocity, and acceleration, identical to the CAA method.Figure 7
Convergence rate of CQ-2x assessed in terms of the maximum relative displacement error.Figure 8
Convergence rate of CQ-2x assessed in terms of the average relative displacement error.Figure 9
Convergence rate of CQ-2x assessed in terms of the maximum relative velocity error.Figure 10
Convergence rate of CQ-2x assessed in terms of the average relative velocity error.Figure 11
Convergence rate of CQ-2x assessed in terms of the maximum relative acceleration error.Figure 12
Convergence rate of CQ-2x assessed in terms of the average relative acceleration error.


## 6.2. Example 2

The purpose of this example is to validate the controllable numerical dissipation and unconditional stability in case of stiffness hardening of the CQ-2x method. In this example, the well-known unforced and undamped Duffing oscillator [28] is considered. The nonlinear dynamics equation is(66)a+S1u1+S2u2=0,where S1=100 and S2−10. The initial conditions are u0=1.5 and v0=0.Equation (66) represents a hardening elastic spring. This is an example of a conservative system that maintains a constant total energy E in its exact solutions. For equation (66),(67)E=12v2+12S1u2+14S1S2u4.To assess the time integration scheme, the percentage error in terms of the energy is introduced [29]:(68)Er=E−E0E0×100%,where E0 is the total energy at t=0.Equation (66) has an exact solution expressed in terms of Jacobian elliptic functions [30]:(69)ut=u0cn−w^t,k,where(70)ω^2=S11+S2u02,k2=S2u0221+S2u02.The periodT of this oscillator is equal to 0.15. The exact phase portrait is given in Figure 13. Also, shown in Figures 14–22 are the phase portraits of numerical results obtained using different integration schemes with the same time step size of Δt=T/20 over a time duration of 100 T. The maximum values of the percentage errors of the nine considered algorithms over this time duration of 100 T are summarized in Table 1 for various time step sizes. When Δt=1/25T, CQ-2xρ∞=1 has an error of 6.06%, whereas the KR-α method has an error of at least 99.63%. It can be seen from Figures 14–17 that ρ∞ is smaller, the numerical damping is bigger, and the undamped Duffing oscillator decays faster. To validate the unconditional stability of CQ-2x in the stiffness hardening case, a time step size of Δt=1/2T over a time duration of 100 T is considered. Figures 23–31 show the phase portraits of the numerical results obtained using the nine integration schemes with the same time step size of Δt=1/2T over a time duration of 100 T. It is meaningless to study an algorithm’s accuracy when Δt=1/2T, but it is useful to explore its stability. Figures 23–26 show that even when Δt=1/2T, CQ-2xρ∞=1,0.8,0.5,0 is still stable. By contrast, the KR-α method is robustly stable only when ρ∞=0, as seen in Figure 30; when ρ∞=1,0.8,and0.5, the KR-α method becomes evidently unstable, as seen in Figures 27–29. e CAA method also exhibits some instability, as seen in Figure 31 (for the CAA method, the iteration termination criterion is that the absolute value of the error between the displacements in the previous and current iterations is less than 10−6).Figure 13
Exact solution for the phase portrait of the hardening spring.Figure 14
Phase portrait of the hardening spring calculated using CQ-2xρ∞=1 when Δt=T/20.Figure 15
Phase portrait of the hardening spring calculated using CQ-2xρ∞=0.8 when Δt=T/20.Figure 16
Phase portrait of the hardening spring calculated using CQ-2xρ∞=0.5 when Δt=T/20.Figure 17
Phase portrait of the hardening spring calculated using CQ-2xρ∞=0 when Δt=T/20.Figure 18
Phase portrait of the hardening spring calculated using the KR-α method ρ∞=1 when Δt=T/20.Figure 19
Phase portrait of the hardening spring calculated using the KR-α method ρ∞=0.8 when Δt=T/20.Figure 20
Phase portrait of the hardening spring calculated using the KR-α method ρ∞=0.5 when Δt=T/20.Figure 21
Phase portrait of the hardening spring calculated using the KR-α method ρ∞=0 when Δt=T/20.Figure 22
Phase portrait of the hardening spring calculated using the CAA method whenΔt=T/20.Table 1
Maximum percentage errors of the total energy over a time duration of 100T for the hardening spring.
CQ-2xKR-αCAAΔtρ∞=1ρ∞=0.8ρ∞=0.5ρ∞=0ρ∞=1ρ∞=0.8ρ∞=0.5ρ∞=0T/10000.0019.8845.5477.2271.6368.5560.601.310.00T/2000.0960.8087.9398.2496.0095.2292.976.830.04T/1000.3679.7696.0799.6899.3298.2597.3615.890.17T/501.4792.1399.0899.97158.3999.3599.0437.140.66T/256.0697.7399.87100.00105670.0299.7499.63400.592.62T/209.5198.6099.94100.00∞99.8199.73700089.874.07T/1516.7799.2899.98100.00∞99.8799.81219047.467.13T/2FiniteFiniteFiniteFinite∞∞∞finite∞Figure 23
Phase portrait of the hardening spring calculated using CQ-2xρ∞=1 when Δt=1/2T.Figure 24
Phase portrait of the hardening spring calculated using CQ-2xρ∞=0.8 when Δt=1/2T.Figure 25
Phase portrait of the hardening spring calculated using CQ-2xρ∞=0.5 when Δt=1/2T.Figure 26
Phase portrait of the hardening spring calculated using CQ-2xρ∞=0 when Δt=1/2T.Figure 27
Phase portrait of the hardening spring calculated using the KR-α method ρ∞=1 when Δt=1/2T.Figure 28
Phase portrait of the hardening spring calculated using the KR-methodρ∞=0.8 when Δt=1/2T.Figure 29
Phase portrait of the hardening spring calculated using the KR-α method ρ∞=0.5 when Δt=1/2T.Figure 30
Phase portrait of the hardening spring calculated using the KR-α method ρ∞=0 when Δt=1/2T.Figure 31
Phase portrait of the hardening spring calculated using the CAA method whenΔt=1/2T.The time consumptions of the different algorithms with different time steps over the time duration of 100T are listed in Table 2. For the same time step size, the noniterative algorithms (CQ-2x and KR-α) are more efficient than the implicit algorithm (CAA). In addition, the CQ-2x method is more efficient than the KR-α method for the same time step size. For example, when Δt=1/25T, the time consumption of CQ-2x is at most 0.0015 seconds, whereas the KR-α method requires at least 0.0048 seconds.Table 2
Time consumptions of the different algorithms for solving the Duffing vibration problem in Example 2 over a time duration of 100T (unit: second).
CQ-2xKR-αCAAΔtρ∞=1ρ∞=0.8ρ∞=0.5ρ∞=0ρ∞=1ρ∞=0.8ρ∞=0.5ρ∞=0T/10000.05900.00900.01800.01740.33050.18710.19360.19471.1631T/2000.01150.00180.00370.00360.06860.04050.04740.04310.2700T/1000.00600.00090.00200.00190.03310.01990.01980.02410.1386T/500.00320.00040.00090.00090.01660.00950.01040.01250.0835T/250.00150.00020.00050.00050.00790.00480.00480.00800.0511T/200.00150.00020.00040.00040.00660.00400.00390.00680.0426T/150.00090.00010.00030.00030.00580.00270.00290.00510.0386


## 6.3. Example 3

The purpose of this example is to validate the efficiency of the CQ-2x method compared with the implicit method and KR−α method. In this example, the forced vibration response of a n-DOF (n = 200, 500, 1000, 2000, and 4000) spring-mass system was considered, as shown in Figure 32. The structural properties of this system are assumed to be mi=100kg and ki=1071−ui−ui−12N/m, in which i=1,2,…,n. This system is subjected to a ground acceleration of 10sinπt. The highest natural frequency is 632.4 rad/s [31]. According to the highest natural frequency, a time step duration of Δt=0.001s is obtained from Newmark’s average acceleration method and is selected to be the exact solution. The responses of displacement and velocity versus time of the 200th DOF when N = 200 are plotted in Figures 33 and 34 . The analysis is accomplished by considering the time step duration equal to 0.02 s. The maximum displacement and velocity error ratios of the 200th DOF when n = 200 are listed in Table 3. It can be seen that the accuracy of the CQ-2x method is sufficient. The maximum error ratio here is defined by the function of max (abs(calculated value-exact value)))/max(abs(exact value)), where max = maximum and abs = absolute. Importantly, without the need for iteration calculations, CQ-2x algorithms have approximately less than half the computational time consumption compared to that consumed by the implicit algorithm CAA. With solving the inverse matrix only once for each time step, CQ-2x algorithms are more efficient than the noniterative KR−α method in which solving the inverse matrix twice for each time step is needed, as listed in Table 4.Figure 32
Schematic plot of then-DOF spring-mass system.Figure 33
Displacement responses of the 200th DOF whenn  = 200.Figure 34
Velocity responses of the 200th DOF whenn  = 200.Table 3
Maximum displacement and velocity responses’ errors of the 200th DOF whenn = 200.
CAACQ-2xKR-αρ∞=1ρ∞=0.8ρ∞=1ρ∞=0.8Displacement0.003310.002840.031170.002830.00289Velocity0.006030.005530.034540.029550.03537Table 4
Time consumptions of the different algorithms for solving then-DOF system in Example 3 (unit: second).
CAACQ-2xKR-αρ∞=1ρ∞=0.8ρ∞=1ρ∞=0.8n=2005.374781.633621.516273.752333.54338n=50039.1792314.8079711.1223131.6716823.96452n=1000246.3178986.0823991.89233182.75176183.34749n=20002045.64400318.19379492.85635893.826941094.05923n=400014992.5984811.7546220.65112205.04014868.966


## 6.4. Example 4

The purpose of this example is to validate the efficiency of the CQ-2x method in case of Vehicle-Bridge Coupled Vibration where the contact nonlinearity exists. In this example, a simple beam of span length = 25 m is subjected to a moving spring mass [32] (as shown in Figure 35). The following data are adopted: Young’s modulus E = 2.87 Gpa, ρ=0.2, moment of inertia I=2.90m4, mass per unit length m=2303kg/m, suspended mass Mv=5750kg, suspension stiffness kv=1595kN/m, and speed v=90km/h. The frequency of vibration computed of the spring mass is ωv=16.66rad/s. The frequencies of vibration computed of the bridge are ω1=30.02rad/s,ω2=120.08rad/s,ω3=270.19rad/s,ω4=480.34rad/s,ω5=750.58rad/s,andω6=1081.0rad/s. The maximum frequency of vibration computed of the bridge is ωmax=95432rad/s. According to stability requirement of the CDM, the widest time step satisfying CDM’s stability is 2/ωmax=2.0957×10−5s. A time step duration of Δt=2×10−5s is obtained for the CDM, and the calculation result is selected to be the exact solution, while with unconditional stability, the time step duration of CQ2xρ∞=1, KR-αρ∞=1, and CAA is only determined by the requirements of the calculation accuracy; here, Δt=0.001s. The midpoint vertical deflection of the beam and deflection of the spring mass versus time are plotted in Figures 36 and 37. The maximum displacement error ratios (defined as that in Example 3) of the midpoint vertical deflection of the beam and the deflection of the spring mass are listed in Table 5. It can be seen that the accuracy of the CQ-2x method is sufficient. The efficiency of the CDM Δt=0.00002s, CQ2xρ∞=1, KR-αρ∞=1, and CAA Δt=0.001s is compared. The average time consumption of the four algorithms in which every algorithm had been performed five times is listed in Table 6. In this example, limited by its conditional stability, the CDM has to be performed with very small Δt=0.00002s, which leads to its efficiency much lower than the other three algorithms. Without the need for iteration calculations, CQ-2x and KR-α algorithms show higher efficiency than that of the implicit algorithm CAA although the maximum number of iteration in each time step of CAA is only 3. The CQ-2x algorithm are more efficient than the noniterative KR−α algorithm since the computation amount in each time step of CQ-2x is much less than that of KR-α.Figure 35
Beam with the moving spring mass.Figure 36
Deflection of the spring mass.Figure 37
Midpoint vertical deflection of the beam.Table 5
Maximum deflection error ratios in Example 4.
CAACQ-2xρ∞=1KR-αρ∞=1Midpoint of beam0.06720.01770.0061Spring mass0.03540.03030.0051Table 6
Time consumption of the different algorithms for solving the beam with the moving spring mass in Example 4 (unit: second).
CDMCAACQ-2xρ∞=1KR-αρ∞=11153.891.409380.512521.21876


## 7. Conclusions

A novel idea based on dimensional analysis for designing noniterative displacement algorithms to be applied in the analysis of nonlinear structural dynamics problems has been introduced and proven to be effective. Originating from only an objective function, a one-parameter family of two-step unconditionally stable noniterative displacement algorithms with controllable numerical energy dissipation has been developed. Analysis results show that the proposed CQ-2x method has unconditional stability in the cases of both stiffness hardening and softening. With the advantages of noniterative formulations of both displacement and velocity, controllable numerical dissipation, the use of one solver per time step, and no overshoot in either displacement or velocity, CQ-2x is expected to serve as a competitive displacement method for solving nonlinear structural dynamics problems.



---

*Source: 4689090-2021-04-21.xml*